From 1f56ed64bda6838ec9fd0988b77d5c71d2d25600 Mon Sep 17 00:00:00 2001
From: andre_ramos <andreramos@lampsco.com>
Date: Fri, 14 Feb 2025 17:01:04 -0300
Subject: [PATCH 1/4] add cycle component

---
 .gitignore                      |   4 +-
 Project.toml                    |   4 +-
 paper_tests/m4_test/m4_test.jl  |   6 +-
 paper_tests/m4_test/m4_test.py  |   7 +-
 src/StateSpaceLearning.jl       |   2 +-
 src/fit_forecast.jl             |  10 +-
 src/models/structural_model.jl  | 356 ++++++++++++++++++++++++++++----
 src/utils.jl                    |  30 ---
 test/models/structural_model.jl | 308 +++++++++++++++++++++++----
 9 files changed, 603 insertions(+), 124 deletions(-)

diff --git a/.gitignore b/.gitignore
index 3e65f9d..151d291 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,4 +13,6 @@ paper_tests/mv_probabilistic_forecast/
 docs/build
 longhorizon/
 test.jl
-test_ms.jl
\ No newline at end of file
+test_ms.jl
+new/
+test2.jl
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index f0a6902..18b24fd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,13 +1,14 @@
 name = "StateSpaceLearning"
 uuid = "971c4b7c-2c4e-4bac-8525-e842df3cde7b"
 authors = ["andreramosfc <andreramosfdc@gmail.com>"]
-version = "1.3.2"
+version = "1.4.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLMNet = "8d5ece8b-de18-5317-b113-243142960cc6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
@@ -20,4 +21,5 @@ PlotsExt = "Plots"
 [compat]
 Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
 GLMNet = "0.5.0, 0.5.1, 0.5.2, 0.6.0, 0.6.1, 0.6.2, 0.6.3, 0.7.0, 0.7.1, 0.7.2, 0.7.3, 0.7.4"
+SparseArrays = "^1.0.0"
 julia = "1"
diff --git a/paper_tests/m4_test/m4_test.jl b/paper_tests/m4_test/m4_test.jl
index 890ed14..0fa958e 100644
--- a/paper_tests/m4_test/m4_test.jl
+++ b/paper_tests/m4_test/m4_test.jl
@@ -25,6 +25,8 @@ function append_results(filepath, results_df)
     if isfile(filepath)
         df_old = CSV.read(filepath, DataFrame)
         results_df = vcat(df_old, results_df)
+        @info "MASE avg = $(mean(results_df[:, :MASE]))"
+        @info "sMAPE avg = $(mean(results_df[:, :sMAPE]))"
     end
     return CSV.write(filepath, results_df)
 end
@@ -138,6 +140,4 @@ end
 
 create_dirs()
 
-main()
-
-#run_config(DataFrame(), false, "aic", 0.1, true, 2794)#max sample size
+main()
\ No newline at end of file
diff --git a/paper_tests/m4_test/m4_test.py b/paper_tests/m4_test/m4_test.py
index cf8836c..d612816 100644
--- a/paper_tests/m4_test/m4_test.py
+++ b/paper_tests/m4_test/m4_test.py
@@ -3,7 +3,12 @@
 import statsmodels.api as sm
 import numpy as np
 
-df_train = pd.read_csv("paper_tests/m4_test/Monthly-train.csv")
+df_train1 = pd.read_csv("paper_tests/m4_test/Monthly-train1.csv")
+df_train2 = pd.read_csv("paper_tests/m4_test/Monthly-train2.csv")
+df_train3 = pd.read_csv("paper_tests/m4_test/Monthly-train3.csv")
+df_train4 = pd.read_csv("paper_tests/m4_test/Monthly-train4.csv")
+df_train = pd.concat([df_train1, df_train2, df_train3, df_train4])
+m4_info = pd.read_csv("paper_tests/m4_test/M4-info.csv")
 df_test = pd.read_csv("paper_tests/m4_test/Monthly-test.csv")
 ssl_init_df = pd.read_csv("paper_tests/m4_test/init_SSL/SSL_aic_0.1_false.csv")
 
diff --git a/src/StateSpaceLearning.jl b/src/StateSpaceLearning.jl
index 047fe6a..d19d07e 100644
--- a/src/StateSpaceLearning.jl
+++ b/src/StateSpaceLearning.jl
@@ -1,6 +1,6 @@
 module StateSpaceLearning
 
-using LinearAlgebra, Statistics, GLMNet, Distributions
+using LinearAlgebra, Statistics, GLMNet, Distributions, SparseArrays
 
 abstract type StateSpaceLearningModel end
 
diff --git a/src/fit_forecast.jl b/src/fit_forecast.jl
index 3b0c7b9..89079a6 100644
--- a/src/fit_forecast.jl
+++ b/src/fit_forecast.jl
@@ -122,15 +122,7 @@ function forecast(
 
     Exogenous_X = model.X[:, exog_idx]
     complete_matrix = create_X(
-        model.level,
-        model.stochastic_level,
-        model.trend,
-        model.stochastic_trend,
-        model.seasonal,
-        model.stochastic_seasonal,
-        model.freq_seasonal,
-        model.outlier,
-        model.ζ_ω_threshold,
+        model,
         Exogenous_X,
         steps_ahead,
         Exogenous_Forecast,
diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl
index c7f301a..3027bf4 100644
--- a/src/models/structural_model.jl
+++ b/src/models/structural_model.jl
@@ -12,6 +12,7 @@ Instantiates a Structural State Space Learning model.
         freq_seasonal::Union{Int, Vector{Int}}=12,
         outlier::Bool=true,
         ζ_ω_threshold::Int=12,
+        stochastic_start::Int=1,
         Exogenous_X::Matrix=if typeof(y) <: Vector
             zeros(length(y), 0)
         else
@@ -59,8 +60,12 @@ mutable struct StructuralModel <: StateSpaceLearningModel
     seasonal::Bool
     stochastic_seasonal::Bool
     freq_seasonal::Union{Int,Vector{Int}}
+    cycle_period::Union{Union{Int, <:AbstractFloat},Vector{Int}, Vector{<:AbstractFloat}}
+    cycle_matrix::Vector{Matrix}
+    stochastic_cycle::Bool
     outlier::Bool
     ζ_ω_threshold::Int
+    stochastic_start::Int
     n_exogenous::Int
     output::Union{Vector{Output},Output,Nothing}
 
@@ -73,8 +78,12 @@ mutable struct StructuralModel <: StateSpaceLearningModel
         seasonal::Bool=true,
         stochastic_seasonal::Bool=true,
         freq_seasonal::Union{Int,Vector{Int}}=12,
+        cycle_period::Union{Union{Int, <:AbstractFloat},Vector{Int}, Vector{<:AbstractFloat}}=0,
+        dumping_cycle::Float64=1.0,
+        stochastic_cycle::Bool=false,
         outlier::Bool=true,
         ζ_ω_threshold::Int=12,
+        stochastic_start::Int=1,
         Exogenous_X::Matrix=if typeof(y) <: Vector
             zeros(length(y), 0)
         else
@@ -88,6 +97,25 @@ mutable struct StructuralModel <: StateSpaceLearningModel
         else
             @assert seasonal ? size(y, 1) > minimum(freq_seasonal) : true "Time series must be longer than the seasonal period"
         end
+        @assert 1 <= stochastic_start < length(y) "stochastic_start must be greater than or equal to 1"
+        @assert 0 < dumping_cycle <= 1 "dumping_cycle must be greater than 0 and less than or equal to 1"
+        if cycle_period != 0 && !isempty(cycle_period)
+            if typeof(cycle_period) <: Vector
+                cycle_matrix = Vector{Matrix}(undef, length(cycle_period))
+                for i in eachindex(cycle_period)
+                    A = dumping_cycle * cos(2 * pi / cycle_period[i])
+                    B = dumping_cycle * sin(2 * pi / cycle_period[i]) 
+                    cycle_matrix[i] = [A B; -B A]
+                end
+            else
+                cycle_matrix = Vector{Matrix}(undef, 1)
+                A = dumping_cycle * cos(2 * pi / cycle_period)
+                B = dumping_cycle * sin(2 * pi / cycle_period) 
+                cycle_matrix[1] = [A B; -B A]
+            end
+        else
+            cycle_matrix = Vector{Matrix}(undef, 0)
+        end
         X = create_X(
             level,
             stochastic_level,
@@ -96,8 +124,11 @@ mutable struct StructuralModel <: StateSpaceLearningModel
             seasonal,
             stochastic_seasonal,
             freq_seasonal,
+            cycle_matrix,
+            stochastic_cycle,
             outlier,
             ζ_ω_threshold,
+            stochastic_start,
             Exogenous_X,
         )
         return new(
@@ -110,8 +141,12 @@ mutable struct StructuralModel <: StateSpaceLearningModel
             seasonal,
             stochastic_seasonal,
             freq_seasonal,
+            cycle_period,
+            cycle_matrix,
+            stochastic_cycle,
             outlier,
             ζ_ω_threshold,
+            stochastic_start,
             n_exogenous,
             nothing,
         )
@@ -119,73 +154,107 @@ mutable struct StructuralModel <: StateSpaceLearningModel
 end
 
 """
-    ξ_size(T::Int)::Int
+    ξ_size(T::Int, stochastic_start::Int)::Int
 
     Calculates the size of ξ innovation matrix based on the input T.
 
     # Arguments
     - `T::Int`: Length of the original time series.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
 
     # Returns
     - `Int`: Size of ξ calculated from T.
 
 """
-ξ_size(T::Int)::Int = T - 2
+ξ_size(T::Int, stochastic_start::Int)::Int = T - max(2, stochastic_start)
 
 """
-ζ_size(T::Int, ζ_ω_threshold::Int)::Int
+ζ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int
 
     Calculates the size of ζ innovation matrix based on the input T.
 
     # Arguments
     - `T::Int`: Length of the original time series.
     - `ζ_ω_threshold::Int`: Stabilize parameter ζ.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
 
     # Returns
     - `Int`: Size of ζ calculated from T.
 
 """
-ζ_size(T::Int, ζ_ω_threshold::Int)::Int = T - ζ_ω_threshold - 2
+ζ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int = max(0, T - ζ_ω_threshold - max(2, stochastic_start))
 
 """
-ω_size(T::Int, s::Int)::Int
+ω_size(T::Int, s::Int, stochastic_start::Int)::Int
 
     Calculates the size of ω innovation matrix based on the input T.
 
     # Arguments
     - `T::Int`: Length of the original time series.
     - `s::Int`: Seasonal period.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
 
     # Returns
     - `Int`: Size of ω calculated from T.
 
 """
-ω_size(T::Int, s::Int, ζ_ω_threshold::Int)::Int = T - ζ_ω_threshold - s + 1
+ω_size(T::Int, s::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int = max(0, T - ζ_ω_threshold - s + 1 - max(0, max(2, stochastic_start) - s))
+
+
+"""
+o_size(T::Int, stochastic_start::Int)::Int
+
+    Calculates the size of outlier matrix based on the input T.
+
+    # Arguments
+    - `T::Int`: Length of the original time series.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
+
+    # Returns
+    - `Int`: Size of o calculated from T.
+
+"""
+o_size(T::Int, stochastic_start::Int)::Int = T - max(1, stochastic_start) + 1
 
 """
-    create_ξ(T::Int, steps_ahead::Int)::Matrix
+    ϕ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int
+
+    Calculates the size of ϕ innovation matrix based on the input T.
+
+    # Arguments
+    - `T::Int`: Length of the original time series.
+    - `ζ_ω_threshold::Int`: Stabilize parameter ζ.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
+
+    # Returns
+    - `Int`: Size of ϕ calculated from T.
+"""
+ϕ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int) = 2 * (T - max(2, stochastic_start) + 1) - (ζ_ω_threshold * 2)
+
+"""
+    create_ξ(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix
 
     Creates a matrix of innovations ξ based on the input sizes, and the desired steps ahead (this is necessary for the forecast function)
 
     # Arguments
     - `T::Int`: Length of the original time series.
     - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0).
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
 
     # Returns
     - `Matrix`: Matrix of innovations ξ constructed based on the input sizes.
 
 """
-function create_ξ(T::Int, steps_ahead::Int)::Matrix
-    ξ_matrix = Matrix{AbstractFloat}(undef, T + steps_ahead, T - 1)
-    for t in 1:(T + steps_ahead)
-        ξ_matrix[t, :] = t < T ? vcat(ones(t - 1), zeros(T - t)) : ones(T - 1)
-    end
-
+function create_ξ(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix
+    stochastic_start = max(2, stochastic_start)
+    ξ_matrix = zeros(T + steps_ahead, T - stochastic_start + 1)
+    ones_indexes = findall(I -> Tuple(I)[1] - (stochastic_start - 2) > Tuple(I)[2], CartesianIndices((T + steps_ahead, T - stochastic_start)))
+    ξ_matrix[ones_indexes] .= 1
     return ξ_matrix[:, 1:(end - 1)]
 end
 
 """
-create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int)::Matrix
+create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
 
     Creates a matrix of innovations ζ based on the input sizes, and the desired steps ahead (this is necessary for the forecast function).
 
@@ -193,25 +262,29 @@ create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int)::Matrix
     - `T::Int`: Length of the original time series.
     - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0).
     - `ζ_ω_threshold::Int`: Stabilize parameter ζ.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
 
     # Returns
     - `Matrix`: Matrix of innovations ζ constructed based on the input sizes.
 
 """
-function create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int)::Matrix
-    ζ_matrix = Matrix{AbstractFloat}(undef, T + steps_ahead, T - 2)
-    for t in 1:(T + steps_ahead)
-        ζ_matrix[t, :] = if t < T
-            vcat(collect((t - 2):-1:1), zeros(T - 2 - length(collect((t - 2):-1:1))))
+function create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
+    stochastic_start = max(2, stochastic_start)
+    ζ_matrix = zeros(T + steps_ahead, T - stochastic_start)
+    
+    for t in 2:(T + steps_ahead)
+        if t < T
+            len = t - stochastic_start
+            ζ_matrix[t, 1:len] .= len:-1:1
         else
-            collect((t - 2):-1:(t - T + 1))
+            ζ_matrix[t, :] .= (t - stochastic_start):-1:(t - T + 1)
         end
     end
     return ζ_matrix[:, 1:(end - ζ_ω_threshold)]
 end
 
 """
-create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int)::Matrix
+create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
 
     Creates a matrix of innovations ω based on the input sizes, and the desired steps ahead (this is necessary for the forecast function).
 
@@ -219,14 +292,18 @@ create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int)::M
     - `T::Int`: Length of the original time series.
     - `freq_seasonal::Int`: Seasonal period.
     - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0).
+    - `ζ_ω_threshold::Int`: Stabilize parameter ζ.
     
     # Returns
     - `Matrix`: Matrix of innovations ω constructed based on the input sizes.
 
 """
-function create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int)::Matrix
-    ω_matrix_size = ω_size(T, freq_seasonal, ζ_ω_threshold) + ζ_ω_threshold
-    ω_matrix = zeros(T + steps_ahead, ω_matrix_size)
+function create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
+    stochastic_start = max(2, stochastic_start)
+    ω_matrix_size = T - freq_seasonal + 1
+    stochastic_start_diff = max(0, stochastic_start - freq_seasonal)
+    
+    ω_matrix = zeros(T + steps_ahead, ω_matrix_size - stochastic_start_diff)
     for t in (freq_seasonal + 1):(T + steps_ahead)
         ωₜ_coefs = zeros(ω_matrix_size)
         Mₜ = Int(floor((t - 1) / freq_seasonal))
@@ -236,11 +313,93 @@ function create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold
             1
         ωₜ_coefs[lag₂[0 .< lag₂ .<= ω_matrix_size + (freq_seasonal - 1)] .- (freq_seasonal - 1)] .=
             -1
-        ω_matrix[t, :] = ωₜ_coefs
+        ω_matrix[t, :] = ωₜ_coefs[1+stochastic_start_diff:end]
     end
     return ω_matrix[:, 1:(end - ζ_ω_threshold)]
 end
 
+
+"""
+create_o_matrix(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix
+
+    Creates a matrix of outliers based on the input sizes, and the desired steps ahead (this is necessary for the forecast function).
+
+    # Arguments
+    - `T::Int`: Length of the original time series.
+    - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0).
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
+    
+    # Returns
+    - `Matrix`: Matrix of outliers constructed based on the input sizes.
+
+"""
+function create_o_matrix(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix
+    stochastic_start = max(1, stochastic_start)
+    rows = stochastic_start:T
+    cols = 1:(T - stochastic_start + 1)
+    values = ones(length(rows))
+    o_matrix = sparse(rows, cols, values, T, T - stochastic_start + 1)
+
+    return vcat(o_matrix, zeros(steps_ahead, length(cols)))
+end
+
+"""
+create_ϕ(X_cycle::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
+
+    Creates a matrix of innovations ϕ based on the input sizes, and the desired steps ahead (this is necessary for the forecast function).
+
+    # Arguments
+    - `X_cycle::Matrix`: deterministic Cycle matrix.
+    - `T::Int`: Length of the original time series.
+    - `steps_ahead::Int64`: Number of steps ahead (for estimation purposes this should be set at 0).
+    - `ζ_ω_threshold::Int`: Stabilize parameter ζ.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
+
+    # Returns
+    - `Matrix`: Matrix of innovations ϕ constructed based on the input sizes.
+"""
+function create_ϕ(c_matrix::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
+
+    num_cols = 2 * (T - stochastic_start + 1)
+    X = Matrix{Float64}(undef, T + steps_ahead, num_cols)
+
+    for (idx, t) in enumerate(stochastic_start:T)
+        X[:, 2 * (idx - 1) + 1] = vcat(zeros(t - 1), c_matrix[1:T - t + 1 + steps_ahead, 1])
+        X[:, 2 * (idx - 1) + 2] = vcat(zeros(t - 1), c_matrix[1:T - t + 1 + steps_ahead, 2])
+    end
+
+    return X[:, 1:end - (ζ_ω_threshold * 2)]
+    
+end
+
+"""
+    create_deterministic_cycle_matrix(cycle_matrix::Vector{Matrix}, T::Int, steps_ahead::Int)::Vector{Matrix}
+
+    Creates a deterministic cycle matrix based on the input parameters.
+
+    # Arguments
+    - `cycle_matrix::Vector{Matrix}`: Vector of cycle matrices.
+    - `T::Int`: Length of the original time series.
+    - `steps_ahead::Int`: Number of steps ahead.
+
+    # Returns
+    - `Vector{Matrix}`: Deterministic cycle matrix constructed based on the input parameters.
+"""
+function create_deterministic_cycle_matrix(cycle_matrix::Vector{Matrix}, T::Int, steps_ahead::Int)::Vector{Matrix}
+    deterministic_cycle_matrix = Vector{Matrix}(undef, length(cycle_matrix))
+    for (idx, c_matrix) in enumerate(cycle_matrix)
+        X_cycle = Matrix{Float64}(undef, T+steps_ahead, 2)
+        cycle_matrix_term = c_matrix ^ 0
+        X_cycle[1, :] = cycle_matrix_term[1, :]
+        for t in 2:T+steps_ahead
+            cycle_matrix_term *= c_matrix
+            X_cycle[t, :] = cycle_matrix_term[1, :]
+        end
+        deterministic_cycle_matrix[idx] = X_cycle
+    end
+    return deterministic_cycle_matrix
+end
+
 """
     create_initial_states_Matrix(
     T::Int, freq_seasonal::Union{Int, Vector{Int}}, steps_ahead::Int, level::Bool, trend::Bool, seasonal::Bool
@@ -267,6 +426,7 @@ function create_initial_states_Matrix(
     level::Bool,
     trend::Bool,
     seasonal::Bool,
+    deterministic_cycle_matrix::Vector{Matrix}
 )::Matrix
     initial_states_matrix = zeros(T + steps_ahead, 0)
     if level
@@ -292,6 +452,12 @@ function create_initial_states_Matrix(
         end
     end
 
+    if !isempty(deterministic_cycle_matrix)
+        for c_matrix in deterministic_cycle_matrix
+            initial_states_matrix = hcat(initial_states_matrix, c_matrix)
+        end
+    end
+
     return initial_states_matrix
 end
 
@@ -306,6 +472,7 @@ create_X(
     freq_seasonal::Union{Int, Vector{Int}},
     outlier::Bool,
     ζ_ω_threshold::Int,
+    stochastic_start::Int
     Exogenous_X::Matrix{Fl},
     steps_ahead::Int=0,
     Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)),
@@ -323,6 +490,7 @@ create_X(
     - `freq_seasonal::Union{Int, Vector{Int}}`: Seasonal period.
     - `outlier::Bool`: Flag for considering outlier component.
     - `ζ_ω_threshold::Int`: Stabilize parameter ζ.
+    - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts.
     - `Exogenous_X::Matrix{Fl}`: Exogenous variables matrix.
     - `steps_ahead::Int`: Number of steps ahead (default: 0).
     - `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast matrix (default: zeros).
@@ -338,17 +506,20 @@ function create_X(
     seasonal::Bool,
     stochastic_seasonal::Bool,
     freq_seasonal::Union{Int,Vector{Int}},
+    cycle_matrix::Vector{Matrix},
+    stochastic_cycle::Bool,
     outlier::Bool,
     ζ_ω_threshold::Int,
+    stochastic_start::Int,
     Exogenous_X::Matrix{Fl},
     steps_ahead::Int=0,
     Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)),
 ) where {Fl<:AbstractFloat,Tl<:AbstractFloat}
     T = size(Exogenous_X, 1)
 
-    ξ_matrix = stochastic_level ? create_ξ(T, steps_ahead) : zeros(T + steps_ahead, 0)
+    ξ_matrix = stochastic_level ? create_ξ(T, steps_ahead, stochastic_start) : zeros(T + steps_ahead, 0)
     ζ_matrix = if stochastic_trend
-        create_ζ(T, steps_ahead, ζ_ω_threshold)
+        create_ζ(T, steps_ahead, ζ_ω_threshold, stochastic_start)
     else
         zeros(T + steps_ahead, 0)
     end
@@ -356,24 +527,79 @@ function create_X(
     ω_matrix = zeros(T + steps_ahead, 0)
     if stochastic_seasonal
         for s in freq_seasonal
-            ω_matrix = hcat(ω_matrix, create_ω(T, s, steps_ahead, ζ_ω_threshold))
+            ω_matrix = hcat(ω_matrix, create_ω(T, s, steps_ahead, ζ_ω_threshold, stochastic_start))
         end
     end
-    o_matrix = outlier ? create_o_matrix(T, steps_ahead) : zeros(T + steps_ahead, 0)
+
+    deterministic_cycle_matrix = create_deterministic_cycle_matrix(cycle_matrix, T, steps_ahead)
+    ϕ_matrix = zeros(T + steps_ahead, 0)
+    if stochastic_cycle
+        for c_matrix in deterministic_cycle_matrix
+            ϕ_matrix = hcat(ϕ_matrix, create_ϕ(c_matrix, T, steps_ahead, ζ_ω_threshold, stochastic_start))
+        end
+    end
+
+    o_matrix = outlier ? create_o_matrix(T, steps_ahead, stochastic_start) : zeros(T + steps_ahead, 0)
 
     initial_states_matrix = create_initial_states_Matrix(
-        T, freq_seasonal, steps_ahead, level, trend, seasonal
+        T, freq_seasonal, steps_ahead, level, trend, seasonal, deterministic_cycle_matrix
     )
     return hcat(
         initial_states_matrix,
         ξ_matrix,
         ζ_matrix,
         ω_matrix,
+        ϕ_matrix,
         o_matrix,
         vcat(Exogenous_X, Exogenous_Forecast),
     )
 end
 
+"""
+create_X(
+    model::StructuralModel,
+    Exogenous_X::Matrix{Fl},
+    steps_ahead::Int=0,
+    Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)),
+) where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+
+    Creates the StateSpaceLearning matrix X based on the model and input parameters.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `Exogenous_X::Matrix{Fl}`: Exogenous variables matrix.
+    - `steps_ahead::Int`: Number of steps ahead (default: 0).
+    - `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast matrix (default: zeros).
+
+    # Returns
+    - `Matrix`: StateSpaceLearning matrix X constructed based on the input parameters.
+"""
+function create_X(
+    model::StructuralModel,
+    Exogenous_X::Matrix{Fl},
+    steps_ahead::Int=0,
+    Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)),
+) where {Fl<:AbstractFloat,Tl<:AbstractFloat}
+    
+    return create_X(
+        model.level,
+        model.stochastic_level,
+        model.trend,
+        model.stochastic_trend,
+        model.seasonal,
+        model.stochastic_seasonal,
+        model.freq_seasonal,
+        model.cycle_matrix,
+        model.stochastic_cycle,
+        model.outlier,
+        model.ζ_ω_threshold,
+        model.stochastic_start,
+        Exogenous_X,
+        steps_ahead,
+        Exogenous_Forecast,
+    )
+end
+
 """
 function get_components_indexes(model::StructuralModel)::Dict
 
@@ -401,7 +627,7 @@ function get_components_indexes(model::StructuralModel)::Dict
     end
 
     if model.trend
-        ν1_indexes = [2]
+        ν1_indexes = [FINAL_INDEX + 1]
         initial_states_indexes = vcat(initial_states_indexes, ν1_indexes)
         FINAL_INDEX += length(ν1_indexes)
     else
@@ -416,12 +642,20 @@ function get_components_indexes(model::StructuralModel)::Dict
             FINAL_INDEX += length(γ_s_indexes)
             push!(γ_indexes, γ_s_indexes)
         end
-    else
-        push!(γ_indexes, Int[])
+    end
+
+    c_indexes = Vector{Int}[]
+    if !isempty(model.cycle_matrix)
+        for _ in eachindex(model.cycle_matrix)
+            c_i_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + 2))
+            initial_states_indexes = vcat(initial_states_indexes, c_i_indexes)
+            FINAL_INDEX += length(c_i_indexes)
+            push!(c_indexes, c_i_indexes)
+        end
     end
 
     if model.stochastic_level
-        ξ_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + ξ_size(T)))
+        ξ_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + ξ_size(T, model.stochastic_start)))
         FINAL_INDEX += length(ξ_indexes)
     else
         ξ_indexes = Int[]
@@ -429,7 +663,7 @@ function get_components_indexes(model::StructuralModel)::Dict
 
     if model.stochastic_trend
         ζ_indexes = collect(
-            (FINAL_INDEX + 1):(FINAL_INDEX + ζ_size(T, model.ζ_ω_threshold))
+            (FINAL_INDEX + 1):(FINAL_INDEX + ζ_size(T, model.ζ_ω_threshold, model.stochastic_start))
         )
         FINAL_INDEX += length(ζ_indexes)
     else
@@ -440,7 +674,7 @@ function get_components_indexes(model::StructuralModel)::Dict
     if model.stochastic_seasonal
         for s in model.freq_seasonal
             ω_s_indexes = collect(
-                (FINAL_INDEX + 1):(FINAL_INDEX + ω_size(T, s, model.ζ_ω_threshold))
+                (FINAL_INDEX + 1):(FINAL_INDEX + ω_size(T, s, model.ζ_ω_threshold, model.stochastic_start))
             )
             FINAL_INDEX += length(ω_s_indexes)
             push!(ω_indexes, ω_s_indexes)
@@ -449,8 +683,21 @@ function get_components_indexes(model::StructuralModel)::Dict
         push!(ω_indexes, Int[])
     end
 
+    ϕ_indexes = Vector{Int}[]
+    if model.stochastic_cycle
+        for _ in eachindex(model.cycle_matrix)
+            ϕ_i_indexes = collect(
+                (FINAL_INDEX + 1):(FINAL_INDEX + ϕ_size(T, model.ζ_ω_threshold, model.stochastic_start))
+            )
+            FINAL_INDEX += length(ϕ_i_indexes)
+            push!(ϕ_indexes, ϕ_i_indexes)
+        end
+    else
+        push!(ϕ_indexes, Int[])
+    end
+
     if model.outlier
-        o_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + o_size(T)))
+        o_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + o_size(T, model.stochastic_start)))
         FINAL_INDEX += length(o_indexes)
     else
         o_indexes = Int[]
@@ -469,10 +716,25 @@ function get_components_indexes(model::StructuralModel)::Dict
     )
 
     for (i, s) in enumerate(model.freq_seasonal)
-        components_indexes_dict["γ1_$s"] = γ_indexes[i]
-        components_indexes_dict["ω_$s"] = ω_indexes[i]
+        if model.seasonal
+            components_indexes_dict["γ1_$s"] = γ_indexes[i]
+        end
+
+        if model.stochastic_seasonal
+            components_indexes_dict["ω_$s"] = ω_indexes[i]
+        end
+        
     end
 
+    if !isempty(model.cycle_matrix)
+        for i in eachindex(model.cycle_period)
+            components_indexes_dict["c1_$i"] = c_indexes[i]
+            if model.stochastic_cycle
+                components_indexes_dict["ϕ_$i"] = ϕ_indexes[i]
+            end
+        end
+    end 
+
     return components_indexes_dict
 end
 
@@ -580,6 +842,12 @@ function get_model_innovations(model::StructuralModel)
             push!(model_innovations, "ω_$s")
         end
     end
+
+    if model.stochastic_cycle
+        for i in eachindex(model.cycle_period)
+            push!(model_innovations, "ϕ_$i")
+        end
+    end
     return model_innovations
 end
 
@@ -600,12 +868,16 @@ function get_innovation_simulation_X(
     model::StructuralModel, innovation::String, steps_ahead::Int
 )
     if innovation == "ξ"
-        return create_ξ(length(model.y) + steps_ahead + 1, 0)
+        return create_ξ(length(model.y) + steps_ahead + 1, 0, model.stochastic_start)
     elseif innovation == "ζ"
-        return create_ζ(length(model.y) + steps_ahead + 1, 0, 1)
+        return create_ζ(length(model.y) + steps_ahead + 1, 0, 1, model.stochastic_start)
     elseif occursin("ω_", innovation)
         s = parse(Int, split(innovation, "_")[2])
-        return create_ω(length(model.y) + steps_ahead + 1, s, 0, 1)
+        return create_ω(length(model.y) + steps_ahead + 1, s, 0, 1, model.stochastic_start)
+    elseif occursin("ϕ_", innovation)
+        i = parse(Int, split(innovation, "_")[2])
+        deterministic_cycle_matrix = create_deterministic_cycle_matrix(model.cycle_matrix, length(model.y), steps_ahead)
+        return create_ϕ(deterministic_cycle_matrix[i], length(model.y), steps_ahead, model.ζ_ω_threshold, model.stochastic_start)
     end
 end
 
diff --git a/src/utils.jl b/src/utils.jl
index c8eba57..93a170a 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -162,36 +162,6 @@ function get_fit_and_residuals(
     end
     return ε_vec, fitted_vec
 end
-"""
-o_size(T::Int)::Int
-
-    Calculates the size of outlier matrix based on the input T.
-
-    # Arguments
-    - `T::Int`: Length of the original time series.
-
-    # Returns
-    - `Int`: Size of o calculated from T.
-
-"""
-o_size(T::Int)::Int = T
-
-"""
-create_o_matrix(T::Int, steps_ahead::Int)::Matrix
-
-    Creates a matrix of outliers based on the input sizes, and the desired steps ahead (this is necessary for the forecast function).
-
-    # Arguments
-    - `T::Int`: Length of the original time series.
-    - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0).
-    
-    # Returns
-    - `Matrix`: Matrix of outliers constructed based on the input sizes.
-
-"""
-function create_o_matrix(T::Int, steps_ahead::Int)::Matrix
-    return vcat(Matrix(1.0 * I, T, T), zeros(steps_ahead, T))
-end
 
 """
     handle_missing_values(
diff --git a/test/models/structural_model.jl b/test/models/structural_model.jl
index 382aba7..b6b45f4 100644
--- a/test/models/structural_model.jl
+++ b/test/models/structural_model.jl
@@ -2,27 +2,107 @@
     y1 = rand(100)
 
     model1 = StateSpaceLearning.StructuralModel(y1)
+    model2 = StateSpaceLearning.StructuralModel(y1; freq_seasonal=[3, 10])
+    model3 = StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2])
+    model4 = StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2], dumping_cycle=0.5)
 
     @test typeof(model1) == StateSpaceLearning.StructuralModel
+    @test typeof(model2) == StateSpaceLearning.StructuralModel
+    @test typeof(model3) == StateSpaceLearning.StructuralModel
+    @test typeof(model4) == StateSpaceLearning.StructuralModel
+
+    @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; stochastic_start=0)
+    @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; stochastic_start=1000)
 
     @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; freq_seasonal=1000)
 
+    @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2], dumping_cycle=-1.0)
+    @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2], dumping_cycle=2.0)
+
     exog_error = ones(100, 3)
     @test_throws AssertionError StateSpaceLearning.StructuralModel(
         y1; Exogenous_X=exog_error
     )
 end
 
-@testset "Innovation matrices" begin
-    @test StateSpaceLearning.ξ_size(10) == 8
-    @test StateSpaceLearning.ζ_size(10, 2) == 6
-    @test StateSpaceLearning.ζ_size(10, 0) == 8
-    @test StateSpaceLearning.ω_size(10, 2, 0) == 9
-    @test StateSpaceLearning.ω_size(10, 2, 2) == 7
-    @test StateSpaceLearning.o_size(10) == 10
+@testset "create_deterministic_cycle_matrix" begin
+    cycle_matrix = Vector{Matrix}(undef, 1)
+    A = 1 * cos(2 * pi / 12)
+    B = 1 * sin(2 * pi / 12) 
+    cycle_matrix[1] = [A B; -B A]
+    det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
+    exp_det_matrix1 = [
+        1.0          0.0;
+        0.866025     0.5;
+        0.5          0.866025;
+        2.77556e-16  1.0;
+       -0.5          0.866025;
+    ]
+
+    n, p = size(det_matrix1[1])
+    for i in 1:n 
+        for j in 1:p
+            @test isapprox(det_matrix1[1][i, j], exp_det_matrix1[i, j]; atol=1e-6)
+        end
+    end
+
+    cycle_period = [3, 12.2]
+    cycle_matrix = Vector{Matrix}(undef, length(cycle_period))
+    for i in eachindex(cycle_period)
+        A = 1 * cos(2 * pi / cycle_period[i])
+        B = 1 * sin(2 * pi / cycle_period[i]) 
+        cycle_matrix[i] = [A B; -B A]
+    end
 
-    X_ξ1 = StateSpaceLearning.create_ξ(5, 0)
-    X_ξ2 = StateSpaceLearning.create_ξ(5, 2)
+    ####
+
+    det_matrix2 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
+    exp_det_matrix2 = [[
+        1.0   0.0
+        -0.5   0.866025
+        -0.5  -0.866025
+         1.0  -6.10623e-16
+        -0.5   0.866025
+    ], 
+    [
+        1.0        0.0;
+        0.870285   0.492548;
+        0.514793   0.857315;
+        0.0257479  0.999668;
+       -0.469977   0.882679;
+    ]
+    ]
+
+    n, p = 5, 2
+    for h in eachindex(det_matrix2)
+        for i in 1:n 
+            for j in 1:p
+                @test isapprox(det_matrix2[h][i, j], exp_det_matrix2[h][i, j]; atol=1e-6)
+            end
+        end
+    end
+end
+
+@testset "Innovation matrices" begin
+    @test StateSpaceLearning.ξ_size(10, 1) == 8
+    @test StateSpaceLearning.ζ_size(10, 2, 1) == 6
+    @test StateSpaceLearning.ζ_size(10, 0, 1) == 8
+    @test StateSpaceLearning.ω_size(10, 2, 0, 1) == 9
+    @test StateSpaceLearning.ω_size(10, 2, 2, 1) == 7
+    @test StateSpaceLearning.o_size(10, 1) == 10
+    @test StateSpaceLearning.ϕ_size(10, 0, 1) == 18
+    @test StateSpaceLearning.ϕ_size(10, 2, 1) == 14
+
+    @test StateSpaceLearning.ξ_size(10, 5) == 5
+    @test StateSpaceLearning.ζ_size(10, 2, 5) == 3
+    @test StateSpaceLearning.ζ_size(10, 0, 5) == 5
+    @test StateSpaceLearning.ω_size(10, 2, 0, 5) == 6
+    @test StateSpaceLearning.ω_size(10, 2, 2, 5) == 4
+    @test StateSpaceLearning.o_size(10, 6) == 5
+    @test StateSpaceLearning.ϕ_size(10, 0, 5) == 12
+
+    X_ξ1 = StateSpaceLearning.create_ξ(5, 0, 1)
+    X_ξ2 = StateSpaceLearning.create_ξ(5, 2, 1)
 
     @test X_ξ1 == [
         0.0 0.0 0.0
@@ -42,9 +122,30 @@ end
         1.0 1.0 1.0
     ]
 
-    X_ζ1 = StateSpaceLearning.create_ζ(5, 0, 0)
-    X_ζ2 = StateSpaceLearning.create_ζ(5, 2, 0)
-    X_ζ3 = StateSpaceLearning.create_ζ(5, 2, 2)
+    X_ξ3 = StateSpaceLearning.create_ξ(5, 0, 3)
+    X_ξ4 = StateSpaceLearning.create_ξ(5, 2, 3)
+
+    @test X_ξ3 == [
+        0.0 0.0
+        0.0 0.0
+        1.0 0.0
+        1.0 1.0
+        1.0 1.0
+    ]
+
+    @test X_ξ4 == [
+        0.0 0.0
+        0.0 0.0
+        1.0 0.0
+        1.0 1.0
+        1.0 1.0
+        1.0 1.0
+        1.0 1.0
+    ]
+
+    X_ζ1 = StateSpaceLearning.create_ζ(5, 0, 0, 1)
+    X_ζ2 = StateSpaceLearning.create_ζ(5, 2, 0, 1)
+    X_ζ3 = StateSpaceLearning.create_ζ(5, 2, 2, 1)
 
     @test X_ζ1 == [
         0.0 0.0 0.0
@@ -78,8 +179,25 @@ end
         1,
     )
 
-    X_ω1 = StateSpaceLearning.create_ω(5, 2, 0, 0)
-    X_ω2 = StateSpaceLearning.create_ω(5, 2, 2, 0)
+    X_ζ4 = StateSpaceLearning.create_ζ(6, 2, 2, 3)
+
+    @test X_ζ4 == reshape(
+        [
+            0.0
+            0.0
+            0.0
+            1.0
+            2.0
+            3.0
+            4.0
+            5.0
+        ],
+        8,
+        1,
+    )
+
+    X_ω1 = StateSpaceLearning.create_ω(5, 2, 0, 0, 1)
+    X_ω2 = StateSpaceLearning.create_ω(5, 2, 2, 0, 1)
 
     @test X_ω1 == [
         0.0 0.0 0.0 0.0
@@ -99,16 +217,48 @@ end
         -1.0 1.0 -1.0 1.0
     ]
 
-    X_o1 = StateSpaceLearning.create_o_matrix(3, 0)
-    X_o2 = StateSpaceLearning.create_o_matrix(3, 2)
+    X_ω3 = StateSpaceLearning.create_ω(5, 2, 0, 0, 3)
+
+    @test X_ω3 == [
+        0.0 0.0 0.0
+        0.0 0.0 0.0
+        1.0 0.0 0.0
+        -1.0 1.0 0.0
+        1.0 -1.0 1.0
+    ]
+
+    X_o1 = StateSpaceLearning.create_o_matrix(3, 0, 1)
+    X_o2 = StateSpaceLearning.create_o_matrix(3, 2, 1)
+    X_o3 = StateSpaceLearning.create_o_matrix(3, 0, 2)
 
     @test X_o1 == Matrix(1.0 * I, 3, 3)
     @test X_o2 == vcat(Matrix(1.0 * I, 3, 3), zeros(2, 3))
+    @test X_o3 == Matrix(1.0 * I, 3, 3)[:, 2:3]
+
+    cycle_matrix = Vector{Matrix}(undef, 1)
+    A = 1 * cos(2 * pi / 12)
+    B = 1 * sin(2 * pi / 12) 
+    cycle_matrix[1] = [A B; -B A]
+    det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
+    X_ϕ1 = StateSpaceLearning.create_ϕ(det_matrix1[1], 5, 0, 0, 1)
+
+    @test all(isapprox.(
+        X_ϕ1, 
+        [
+            1.0          0.0       0.0          0.0       0.0       0.0       0.0       0.0  0.0  0.0;
+            0.866025     0.5       1.0          0.0       0.0       0.0       0.0       0.0  0.0  0.0;
+            0.5          0.866025  0.866025     0.5       1.0       0.0       0.0       0.0  0.0  0.0;
+            2.77556e-16  1.0       0.5          0.866025  0.866025  0.5       1.0       0.0  0.0  0.0;
+           -0.5          0.866025  2.77556e-16  1.0       0.5       0.866025  0.866025  0.5  1.0  0.0
+        ],
+        atol=1e-6
+    ))
 end
 
 @testset "Initial State Matrix" begin
-    X1 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, true)
-    X2 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, true)
+    
+    X1 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, true, Vector{Matrix}(undef, 0))
+    X2 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, true, Vector{Matrix}(undef, 0))
 
     @test X1 == [
         1.0 0.0 1.0 0.0
@@ -128,8 +278,8 @@ end
         1.0 6.0 1.0 0.0
     ]
 
-    X3 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, false)
-    X4 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, false)
+    X3 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, false, Vector{Matrix}(undef, 0))
+    X4 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, false, Vector{Matrix}(undef, 0))
 
     @test X3 == [
         1.0 0.0
@@ -149,17 +299,53 @@ end
         1.0 6.0
     ]
 
-    X5 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, false, false)
-    X6 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, false, false)
+    X5 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, false, false, Vector{Matrix}(undef, 0))
+    X6 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, false, false, Vector{Matrix}(undef, 0))
 
     @test X5 == ones(5, 1)
     @test X6 == ones(7, 1)
 
-    X7 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, false, true, false)
-    X8 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, false, true, false)
+    X7 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, false, true, false, Vector{Matrix}(undef, 0))
+    X8 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, false, true, false, Vector{Matrix}(undef, 0))
 
     @test X7 == [0.0; 1.0; 2.0; 3.0; 4.0][:, :]
     @test X8 == [0.0; 1.0; 2.0; 3.0; 4.0; 5.0; 6.0][:, :]
+
+    cycle_matrix = Vector{Matrix}(undef, 1)
+    A = 1 * cos(2 * pi / 12)
+    B = 1 * sin(2 * pi / 12) 
+    cycle_matrix[1] = [A B; -B A]
+    det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
+    det_matrix2 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 2)
+
+    X9 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, true, det_matrix1)
+    X10 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, true, det_matrix2)
+
+    @test all(isapprox.(
+        X9, 
+        [
+            1.0  0.0  1.0  0.0   1.0          0.0;
+            1.0  1.0  0.0  1.0   0.866025     0.5;
+            1.0  2.0  1.0  0.0   0.5          0.866025;
+            1.0  3.0  0.0  1.0   2.77556e-16  1.0;
+            1.0  4.0  1.0  0.0  -0.5          0.866025;
+        ],
+        atol=1e-6
+    ))
+
+    @test all(isapprox.(
+        X10, 
+        [
+            1.0  0.0  1.0  0.0   1.0          0.0;
+            1.0  1.0  0.0  1.0   0.866025     0.5;
+            1.0  2.0  1.0  0.0   0.5          0.866025;
+            1.0  3.0  0.0  1.0   2.77556e-16  1.0;
+            1.0  4.0  1.0  0.0  -0.5          0.866025;
+            1.0  5.0  0.0  1.0  -0.866025     0.5;
+            1.0  6.0  1.0  0.0  -1.0          5.55112e-16
+        ],
+        atol=1e-6
+    ))
 end
 
 @testset "Create X matrix" begin
@@ -168,6 +354,10 @@ end
 
     Exogenous_forecast1 = rand(2, 3)
 
+    cycle_matrix = Vector{Matrix}(undef, 0)
+    stochastic_cycle = false
+    stochastic_start = 1
+
     param_combination = [
         Any[true, 0, 0],
         Any[true, 2, 0],
@@ -231,23 +421,23 @@ end
     ]
     counter = 1
     for args in [
-        [true, true, true, true, true, true, 2, true, 12],
-        [true, true, false, false, false, false, 2, true, 12],
-        [true, true, true, true, false, false, 2, true, 12],
-        [true, false, true, true, false, false, 2, true, 12],
+        [true, true, true, true, true, true, 2, cycle_matrix, stochastic_cycle, true, 12],
+        [true, true, false, false, false, false, 2, cycle_matrix, stochastic_cycle, true, 12],
+        [true, true, true, true, false, false, 2, cycle_matrix, stochastic_cycle, true, 12],
+        [true, false, true, true, false, false, 2, cycle_matrix, stochastic_cycle, true, 12],
     ]
         args = [x in [0, 1] ? Bool(x) : x for x in args]
         for param in param_combination
-            args[8] = param[1]
-            args[9] = param[2]
+            args[end - 1] = param[1]
+            args[end] = param[2]
             if param[3] != 0
                 X1 = StateSpaceLearning.create_X(
-                    args..., Exogenous_X1, param[3], Exogenous_forecast1
+                    args..., stochastic_start, Exogenous_X1, param[3], Exogenous_forecast1
                 )
             else
-                X1 = StateSpaceLearning.create_X(args..., Exogenous_X1, param[3])
+                X1 = StateSpaceLearning.create_X(args..., stochastic_start, Exogenous_X1, param[3])
             end
-            X2 = StateSpaceLearning.create_X(args..., Exogenous_X2, param[3])
+            X2 = StateSpaceLearning.create_X(args..., stochastic_start, Exogenous_X2, param[3])
             @test size(X1) == size_vec1[counter]
             @test size(X2) == size_vec2[counter]
             counter += 1
@@ -391,11 +581,26 @@ end
         ζ_ω_threshold=0,
         Exogenous_X=Exogenous_X2,
     )
+    Local_Linear_Trend_cycle = StateSpaceLearning.StructuralModel(
+        rand(10);
+        level=true,
+        stochastic_level=true,
+        trend=true,
+        stochastic_trend=true,
+        seasonal=false,
+        stochastic_seasonal=false,
+        freq_seasonal=2,
+        cycle_period=3,
+        stochastic_cycle=true,
+        outlier=false,
+        ζ_ω_threshold=0,
+        Exogenous_X=Exogenous_X2,
+    )
 
-    models = [Basic_Structural, Basic_Structural2, Local_Level, Local_Linear_Trend]
+    models = [Basic_Structural, Basic_Structural2, Local_Level, Local_Linear_Trend, Local_Linear_Trend_cycle]
 
     params_vec = [
-        ["ξ", "ζ", "ω_2", "ε"], ["ξ", "ζ", "ω_2", "ω_5", "ε"], ["ξ", "ε"], ["ξ", "ζ", "ε"]
+        ["ξ", "ζ", "ω_2", "ε"], ["ξ", "ζ", "ω_2", "ω_5", "ε"], ["ξ", "ε"], ["ξ", "ζ", "ε"], ["ξ", "ζ", "ε", "ϕ_1"]
     ]
 
     for idx in eachindex(models)
@@ -474,11 +679,26 @@ end
         ζ_ω_threshold=0,
         Exogenous_X=zeros(10, 0),
     )
+    model6 = StateSpaceLearning.StructuralModel(
+        rand(10);
+        level=true,
+        stochastic_level=true,
+        trend=true,
+        stochastic_trend=false,
+        seasonal=true,
+        stochastic_seasonal=true,
+        freq_seasonal=[2, 5],
+        cycle_period=3,
+        stochastic_cycle=true,
+        outlier=true,
+        ζ_ω_threshold=0,
+        Exogenous_X=zeros(10, 0),
+    )
 
     models = [model1, model2, model3, model4]
 
     keys_vec = [
-        ["ξ", "ζ", "ω_2"], ["ζ", "ω_2"], ["ξ", "ω_2"], ["ξ", "ζ"], ["ξ", "ω_2", "ω_5"]
+        ["ξ", "ζ", "ω_2"], ["ζ", "ω_2"], ["ξ", "ω_2"], ["ξ", "ζ"], ["ξ", "ω_2", "ω_5", "ϕ_1"]
     ]
 
     for idx in eachindex(models)
@@ -492,6 +712,7 @@ end
     innovation1 = "ξ"
     innovation2 = "ζ"
     innovation3 = "ω_2"
+    innovation4 = "ϕ_1"
 
     model = StateSpaceLearning.StructuralModel(
         rand(3);
@@ -503,6 +724,8 @@ end
         stochastic_seasonal=true,
         freq_seasonal=2,
         outlier=true,
+        cycle_period=3,
+        stochastic_cycle=true,
         ζ_ω_threshold=0,
         Exogenous_X=zeros(10, 0),
     )
@@ -537,4 +760,17 @@ end
         -1.0 1.0 -1.0 1.0
         0.0 -1.0 1.0 -1.0
     ]
+
+    X4 = StateSpaceLearning.get_innovation_simulation_X(model, innovation4, steps_ahead)
+    @assert all(isapprox.(
+        X4, 
+        [
+            1.0   0.0           0.0   0.0           0.0   0.0;
+            -0.5   0.866025      1.0   0.0           0.0   0.0;
+            -0.5  -0.866025     -0.5   0.866025      1.0   0.0;
+             1.0  -6.10623e-16  -0.5  -0.866025     -0.5   0.866025;
+            -0.5   0.866025      1.0  -6.10623e-16  -0.5  -0.866025;
+        ],
+        atol=1e-6
+    ))
 end

From cd041527fdec496db8f79df18e63b473a869a4d6 Mon Sep 17 00:00:00 2001
From: andre_ramos <andreramos@lampsco.com>
Date: Fri, 14 Feb 2025 17:05:26 -0300
Subject: [PATCH 2/4] format in bluestyle

---
 paper_tests/m4_test/m4_test.jl  |   2 +-
 src/fit_forecast.jl             |   7 +-
 src/models/structural_model.jl  | 127 ++++++++++-----
 test/models/structural_model.jl | 264 +++++++++++++++++++++-----------
 4 files changed, 259 insertions(+), 141 deletions(-)

diff --git a/paper_tests/m4_test/m4_test.jl b/paper_tests/m4_test/m4_test.jl
index 0fa958e..cda5c3a 100644
--- a/paper_tests/m4_test/m4_test.jl
+++ b/paper_tests/m4_test/m4_test.jl
@@ -140,4 +140,4 @@ end
 
 create_dirs()
 
-main()
\ No newline at end of file
+main()
diff --git a/src/fit_forecast.jl b/src/fit_forecast.jl
index 89079a6..5836d7f 100644
--- a/src/fit_forecast.jl
+++ b/src/fit_forecast.jl
@@ -121,12 +121,7 @@ function forecast(
     @assert size(Exogenous_Forecast, 1) == steps_ahead "Exogenous_Forecast must have the same number of rows as steps_ahead"
 
     Exogenous_X = model.X[:, exog_idx]
-    complete_matrix = create_X(
-        model,
-        Exogenous_X,
-        steps_ahead,
-        Exogenous_Forecast,
-    )
+    complete_matrix = create_X(model, Exogenous_X, steps_ahead, Exogenous_Forecast)
 
     if typeof(model.output) == Output
         return AbstractFloat.(
diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl
index 3027bf4..b1d0f5a 100644
--- a/src/models/structural_model.jl
+++ b/src/models/structural_model.jl
@@ -60,7 +60,7 @@ mutable struct StructuralModel <: StateSpaceLearningModel
     seasonal::Bool
     stochastic_seasonal::Bool
     freq_seasonal::Union{Int,Vector{Int}}
-    cycle_period::Union{Union{Int, <:AbstractFloat},Vector{Int}, Vector{<:AbstractFloat}}
+    cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}
     cycle_matrix::Vector{Matrix}
     stochastic_cycle::Bool
     outlier::Bool
@@ -78,7 +78,7 @@ mutable struct StructuralModel <: StateSpaceLearningModel
         seasonal::Bool=true,
         stochastic_seasonal::Bool=true,
         freq_seasonal::Union{Int,Vector{Int}}=12,
-        cycle_period::Union{Union{Int, <:AbstractFloat},Vector{Int}, Vector{<:AbstractFloat}}=0,
+        cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}=0,
         dumping_cycle::Float64=1.0,
         stochastic_cycle::Bool=false,
         outlier::Bool=true,
@@ -104,13 +104,13 @@ mutable struct StructuralModel <: StateSpaceLearningModel
                 cycle_matrix = Vector{Matrix}(undef, length(cycle_period))
                 for i in eachindex(cycle_period)
                     A = dumping_cycle * cos(2 * pi / cycle_period[i])
-                    B = dumping_cycle * sin(2 * pi / cycle_period[i]) 
+                    B = dumping_cycle * sin(2 * pi / cycle_period[i])
                     cycle_matrix[i] = [A B; -B A]
                 end
             else
                 cycle_matrix = Vector{Matrix}(undef, 1)
                 A = dumping_cycle * cos(2 * pi / cycle_period)
-                B = dumping_cycle * sin(2 * pi / cycle_period) 
+                B = dumping_cycle * sin(2 * pi / cycle_period)
                 cycle_matrix[1] = [A B; -B A]
             end
         else
@@ -182,7 +182,8 @@ end
     - `Int`: Size of ζ calculated from T.
 
 """
-ζ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int = max(0, T - ζ_ω_threshold - max(2, stochastic_start))
+ζ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int =
+    max(0, T - ζ_ω_threshold - max(2, stochastic_start))
 
 """
 ω_size(T::Int, s::Int, stochastic_start::Int)::Int
@@ -198,8 +199,8 @@ end
     - `Int`: Size of ω calculated from T.
 
 """
-ω_size(T::Int, s::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int = max(0, T - ζ_ω_threshold - s + 1 - max(0, max(2, stochastic_start) - s))
-
+ω_size(T::Int, s::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int =
+    max(0, T - ζ_ω_threshold - s + 1 - max(0, max(2, stochastic_start) - s))
 
 """
 o_size(T::Int, stochastic_start::Int)::Int
@@ -229,7 +230,8 @@ o_size(T::Int, stochastic_start::Int)::Int = T - max(1, stochastic_start) + 1
     # Returns
     - `Int`: Size of ϕ calculated from T.
 """
-ϕ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int) = 2 * (T - max(2, stochastic_start) + 1) - (ζ_ω_threshold * 2)
+ϕ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int) =
+    2 * (T - max(2, stochastic_start) + 1) - (ζ_ω_threshold * 2)
 
 """
     create_ξ(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix
@@ -248,7 +250,10 @@ o_size(T::Int, stochastic_start::Int)::Int = T - max(1, stochastic_start) + 1
 function create_ξ(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix
     stochastic_start = max(2, stochastic_start)
     ξ_matrix = zeros(T + steps_ahead, T - stochastic_start + 1)
-    ones_indexes = findall(I -> Tuple(I)[1] - (stochastic_start - 2) > Tuple(I)[2], CartesianIndices((T + steps_ahead, T - stochastic_start)))
+    ones_indexes = findall(
+        I -> Tuple(I)[1] - (stochastic_start - 2) > Tuple(I)[2],
+        CartesianIndices((T + steps_ahead, T - stochastic_start)),
+    )
     ξ_matrix[ones_indexes] .= 1
     return ξ_matrix[:, 1:(end - 1)]
 end
@@ -268,10 +273,12 @@ create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)
     - `Matrix`: Matrix of innovations ζ constructed based on the input sizes.
 
 """
-function create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
+function create_ζ(
+    T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int
+)::Matrix
     stochastic_start = max(2, stochastic_start)
     ζ_matrix = zeros(T + steps_ahead, T - stochastic_start)
-    
+
     for t in 2:(T + steps_ahead)
         if t < T
             len = t - stochastic_start
@@ -298,11 +305,13 @@ create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int, st
     - `Matrix`: Matrix of innovations ω constructed based on the input sizes.
 
 """
-function create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
+function create_ω(
+    T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int
+)::Matrix
     stochastic_start = max(2, stochastic_start)
     ω_matrix_size = T - freq_seasonal + 1
     stochastic_start_diff = max(0, stochastic_start - freq_seasonal)
-    
+
     ω_matrix = zeros(T + steps_ahead, ω_matrix_size - stochastic_start_diff)
     for t in (freq_seasonal + 1):(T + steps_ahead)
         ωₜ_coefs = zeros(ω_matrix_size)
@@ -313,12 +322,11 @@ function create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold
             1
         ωₜ_coefs[lag₂[0 .< lag₂ .<= ω_matrix_size + (freq_seasonal - 1)] .- (freq_seasonal - 1)] .=
             -1
-        ω_matrix[t, :] = ωₜ_coefs[1+stochastic_start_diff:end]
+        ω_matrix[t, :] = ωₜ_coefs[(1 + stochastic_start_diff):end]
     end
     return ω_matrix[:, 1:(end - ζ_ω_threshold)]
 end
 
-
 """
 create_o_matrix(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix
 
@@ -358,18 +366,22 @@ create_ϕ(X_cycle::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, sto
     # Returns
     - `Matrix`: Matrix of innovations ϕ constructed based on the input sizes.
 """
-function create_ϕ(c_matrix::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
-
+function create_ϕ(
+    c_matrix::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, stochastic_start::Int
+)::Matrix
     num_cols = 2 * (T - stochastic_start + 1)
     X = Matrix{Float64}(undef, T + steps_ahead, num_cols)
 
     for (idx, t) in enumerate(stochastic_start:T)
-        X[:, 2 * (idx - 1) + 1] = vcat(zeros(t - 1), c_matrix[1:T - t + 1 + steps_ahead, 1])
-        X[:, 2 * (idx - 1) + 2] = vcat(zeros(t - 1), c_matrix[1:T - t + 1 + steps_ahead, 2])
+        X[:, 2 * (idx - 1) + 1] = vcat(
+            zeros(t - 1), c_matrix[1:(T - t + 1 + steps_ahead), 1]
+        )
+        X[:, 2 * (idx - 1) + 2] = vcat(
+            zeros(t - 1), c_matrix[1:(T - t + 1 + steps_ahead), 2]
+        )
     end
 
-    return X[:, 1:end - (ζ_ω_threshold * 2)]
-    
+    return X[:, 1:(end - (ζ_ω_threshold * 2))]
 end
 
 """
@@ -385,13 +397,15 @@ end
     # Returns
     - `Vector{Matrix}`: Deterministic cycle matrix constructed based on the input parameters.
 """
-function create_deterministic_cycle_matrix(cycle_matrix::Vector{Matrix}, T::Int, steps_ahead::Int)::Vector{Matrix}
+function create_deterministic_cycle_matrix(
+    cycle_matrix::Vector{Matrix}, T::Int, steps_ahead::Int
+)::Vector{Matrix}
     deterministic_cycle_matrix = Vector{Matrix}(undef, length(cycle_matrix))
     for (idx, c_matrix) in enumerate(cycle_matrix)
-        X_cycle = Matrix{Float64}(undef, T+steps_ahead, 2)
-        cycle_matrix_term = c_matrix ^ 0
+        X_cycle = Matrix{Float64}(undef, T + steps_ahead, 2)
+        cycle_matrix_term = c_matrix^0
         X_cycle[1, :] = cycle_matrix_term[1, :]
-        for t in 2:T+steps_ahead
+        for t in 2:(T + steps_ahead)
             cycle_matrix_term *= c_matrix
             X_cycle[t, :] = cycle_matrix_term[1, :]
         end
@@ -426,7 +440,7 @@ function create_initial_states_Matrix(
     level::Bool,
     trend::Bool,
     seasonal::Bool,
-    deterministic_cycle_matrix::Vector{Matrix}
+    deterministic_cycle_matrix::Vector{Matrix},
 )::Matrix
     initial_states_matrix = zeros(T + steps_ahead, 0)
     if level
@@ -517,7 +531,11 @@ function create_X(
 ) where {Fl<:AbstractFloat,Tl<:AbstractFloat}
     T = size(Exogenous_X, 1)
 
-    ξ_matrix = stochastic_level ? create_ξ(T, steps_ahead, stochastic_start) : zeros(T + steps_ahead, 0)
+    ξ_matrix = if stochastic_level
+        create_ξ(T, steps_ahead, stochastic_start)
+    else
+        zeros(T + steps_ahead, 0)
+    end
     ζ_matrix = if stochastic_trend
         create_ζ(T, steps_ahead, ζ_ω_threshold, stochastic_start)
     else
@@ -527,19 +545,30 @@ function create_X(
     ω_matrix = zeros(T + steps_ahead, 0)
     if stochastic_seasonal
         for s in freq_seasonal
-            ω_matrix = hcat(ω_matrix, create_ω(T, s, steps_ahead, ζ_ω_threshold, stochastic_start))
+            ω_matrix = hcat(
+                ω_matrix, create_ω(T, s, steps_ahead, ζ_ω_threshold, stochastic_start)
+            )
         end
     end
 
-    deterministic_cycle_matrix = create_deterministic_cycle_matrix(cycle_matrix, T, steps_ahead)
+    deterministic_cycle_matrix = create_deterministic_cycle_matrix(
+        cycle_matrix, T, steps_ahead
+    )
     ϕ_matrix = zeros(T + steps_ahead, 0)
     if stochastic_cycle
         for c_matrix in deterministic_cycle_matrix
-            ϕ_matrix = hcat(ϕ_matrix, create_ϕ(c_matrix, T, steps_ahead, ζ_ω_threshold, stochastic_start))
+            ϕ_matrix = hcat(
+                ϕ_matrix,
+                create_ϕ(c_matrix, T, steps_ahead, ζ_ω_threshold, stochastic_start),
+            )
         end
     end
 
-    o_matrix = outlier ? create_o_matrix(T, steps_ahead, stochastic_start) : zeros(T + steps_ahead, 0)
+    o_matrix = if outlier
+        create_o_matrix(T, steps_ahead, stochastic_start)
+    else
+        zeros(T + steps_ahead, 0)
+    end
 
     initial_states_matrix = create_initial_states_Matrix(
         T, freq_seasonal, steps_ahead, level, trend, seasonal, deterministic_cycle_matrix
@@ -580,7 +609,6 @@ function create_X(
     steps_ahead::Int=0,
     Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)),
 ) where {Fl<:AbstractFloat,Tl<:AbstractFloat}
-    
     return create_X(
         model.level,
         model.stochastic_level,
@@ -655,7 +683,9 @@ function get_components_indexes(model::StructuralModel)::Dict
     end
 
     if model.stochastic_level
-        ξ_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + ξ_size(T, model.stochastic_start)))
+        ξ_indexes = collect(
+            (FINAL_INDEX + 1):(FINAL_INDEX + ξ_size(T, model.stochastic_start))
+        )
         FINAL_INDEX += length(ξ_indexes)
     else
         ξ_indexes = Int[]
@@ -663,7 +693,9 @@ function get_components_indexes(model::StructuralModel)::Dict
 
     if model.stochastic_trend
         ζ_indexes = collect(
-            (FINAL_INDEX + 1):(FINAL_INDEX + ζ_size(T, model.ζ_ω_threshold, model.stochastic_start))
+            (FINAL_INDEX + 1):(FINAL_INDEX + ζ_size(
+                T, model.ζ_ω_threshold, model.stochastic_start
+            )),
         )
         FINAL_INDEX += length(ζ_indexes)
     else
@@ -674,7 +706,9 @@ function get_components_indexes(model::StructuralModel)::Dict
     if model.stochastic_seasonal
         for s in model.freq_seasonal
             ω_s_indexes = collect(
-                (FINAL_INDEX + 1):(FINAL_INDEX + ω_size(T, s, model.ζ_ω_threshold, model.stochastic_start))
+                (FINAL_INDEX + 1):(FINAL_INDEX + ω_size(
+                    T, s, model.ζ_ω_threshold, model.stochastic_start
+                )),
             )
             FINAL_INDEX += length(ω_s_indexes)
             push!(ω_indexes, ω_s_indexes)
@@ -687,7 +721,9 @@ function get_components_indexes(model::StructuralModel)::Dict
     if model.stochastic_cycle
         for _ in eachindex(model.cycle_matrix)
             ϕ_i_indexes = collect(
-                (FINAL_INDEX + 1):(FINAL_INDEX + ϕ_size(T, model.ζ_ω_threshold, model.stochastic_start))
+                (FINAL_INDEX + 1):(FINAL_INDEX + ϕ_size(
+                    T, model.ζ_ω_threshold, model.stochastic_start
+                )),
             )
             FINAL_INDEX += length(ϕ_i_indexes)
             push!(ϕ_indexes, ϕ_i_indexes)
@@ -697,7 +733,9 @@ function get_components_indexes(model::StructuralModel)::Dict
     end
 
     if model.outlier
-        o_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + o_size(T, model.stochastic_start)))
+        o_indexes = collect(
+            (FINAL_INDEX + 1):(FINAL_INDEX + o_size(T, model.stochastic_start))
+        )
         FINAL_INDEX += length(o_indexes)
     else
         o_indexes = Int[]
@@ -723,7 +761,6 @@ function get_components_indexes(model::StructuralModel)::Dict
         if model.stochastic_seasonal
             components_indexes_dict["ω_$s"] = ω_indexes[i]
         end
-        
     end
 
     if !isempty(model.cycle_matrix)
@@ -733,7 +770,7 @@ function get_components_indexes(model::StructuralModel)::Dict
                 components_indexes_dict["ϕ_$i"] = ϕ_indexes[i]
             end
         end
-    end 
+    end
 
     return components_indexes_dict
 end
@@ -876,8 +913,16 @@ function get_innovation_simulation_X(
         return create_ω(length(model.y) + steps_ahead + 1, s, 0, 1, model.stochastic_start)
     elseif occursin("ϕ_", innovation)
         i = parse(Int, split(innovation, "_")[2])
-        deterministic_cycle_matrix = create_deterministic_cycle_matrix(model.cycle_matrix, length(model.y), steps_ahead)
-        return create_ϕ(deterministic_cycle_matrix[i], length(model.y), steps_ahead, model.ζ_ω_threshold, model.stochastic_start)
+        deterministic_cycle_matrix = create_deterministic_cycle_matrix(
+            model.cycle_matrix, length(model.y), steps_ahead
+        )
+        return create_ϕ(
+            deterministic_cycle_matrix[i],
+            length(model.y),
+            steps_ahead,
+            model.ζ_ω_threshold,
+            model.stochastic_start,
+        )
     end
 end
 
diff --git a/test/models/structural_model.jl b/test/models/structural_model.jl
index b6b45f4..fe7a091 100644
--- a/test/models/structural_model.jl
+++ b/test/models/structural_model.jl
@@ -4,7 +4,9 @@
     model1 = StateSpaceLearning.StructuralModel(y1)
     model2 = StateSpaceLearning.StructuralModel(y1; freq_seasonal=[3, 10])
     model3 = StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2])
-    model4 = StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2], dumping_cycle=0.5)
+    model4 = StateSpaceLearning.StructuralModel(
+        y1; cycle_period=[3, 10.2], dumping_cycle=0.5
+    )
 
     @test typeof(model1) == StateSpaceLearning.StructuralModel
     @test typeof(model2) == StateSpaceLearning.StructuralModel
@@ -12,12 +14,18 @@
     @test typeof(model4) == StateSpaceLearning.StructuralModel
 
     @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; stochastic_start=0)
-    @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; stochastic_start=1000)
+    @test_throws AssertionError StateSpaceLearning.StructuralModel(
+        y1; stochastic_start=1000
+    )
 
     @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; freq_seasonal=1000)
 
-    @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2], dumping_cycle=-1.0)
-    @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2], dumping_cycle=2.0)
+    @test_throws AssertionError StateSpaceLearning.StructuralModel(
+        y1; cycle_period=[3, 10.2], dumping_cycle=-1.0
+    )
+    @test_throws AssertionError StateSpaceLearning.StructuralModel(
+        y1; cycle_period=[3, 10.2], dumping_cycle=2.0
+    )
 
     exog_error = ones(100, 3)
     @test_throws AssertionError StateSpaceLearning.StructuralModel(
@@ -28,19 +36,19 @@ end
 @testset "create_deterministic_cycle_matrix" begin
     cycle_matrix = Vector{Matrix}(undef, 1)
     A = 1 * cos(2 * pi / 12)
-    B = 1 * sin(2 * pi / 12) 
+    B = 1 * sin(2 * pi / 12)
     cycle_matrix[1] = [A B; -B A]
     det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
     exp_det_matrix1 = [
-        1.0          0.0;
-        0.866025     0.5;
-        0.5          0.866025;
-        2.77556e-16  1.0;
-       -0.5          0.866025;
+        1.0 0.0
+        0.866025 0.5
+        0.5 0.866025
+        2.77556e-16 1.0
+        -0.5 0.866025
     ]
 
     n, p = size(det_matrix1[1])
-    for i in 1:n 
+    for i in 1:n
         for j in 1:p
             @test isapprox(det_matrix1[1][i, j], exp_det_matrix1[i, j]; atol=1e-6)
         end
@@ -50,32 +58,33 @@ end
     cycle_matrix = Vector{Matrix}(undef, length(cycle_period))
     for i in eachindex(cycle_period)
         A = 1 * cos(2 * pi / cycle_period[i])
-        B = 1 * sin(2 * pi / cycle_period[i]) 
+        B = 1 * sin(2 * pi / cycle_period[i])
         cycle_matrix[i] = [A B; -B A]
     end
 
     ####
 
     det_matrix2 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
-    exp_det_matrix2 = [[
-        1.0   0.0
-        -0.5   0.866025
-        -0.5  -0.866025
-         1.0  -6.10623e-16
-        -0.5   0.866025
-    ], 
-    [
-        1.0        0.0;
-        0.870285   0.492548;
-        0.514793   0.857315;
-        0.0257479  0.999668;
-       -0.469977   0.882679;
-    ]
+    exp_det_matrix2 = [
+        [
+            1.0 0.0
+            -0.5 0.866025
+            -0.5 -0.866025
+            1.0 -6.10623e-16
+            -0.5 0.866025
+        ],
+        [
+            1.0 0.0
+            0.870285 0.492548
+            0.514793 0.857315
+            0.0257479 0.999668
+            -0.469977 0.882679
+        ],
     ]
 
     n, p = 5, 2
     for h in eachindex(det_matrix2)
-        for i in 1:n 
+        for i in 1:n
             for j in 1:p
                 @test isapprox(det_matrix2[h][i, j], exp_det_matrix2[h][i, j]; atol=1e-6)
             end
@@ -237,28 +246,33 @@ end
 
     cycle_matrix = Vector{Matrix}(undef, 1)
     A = 1 * cos(2 * pi / 12)
-    B = 1 * sin(2 * pi / 12) 
+    B = 1 * sin(2 * pi / 12)
     cycle_matrix[1] = [A B; -B A]
     det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
     X_ϕ1 = StateSpaceLearning.create_ϕ(det_matrix1[1], 5, 0, 0, 1)
 
-    @test all(isapprox.(
-        X_ϕ1, 
-        [
-            1.0          0.0       0.0          0.0       0.0       0.0       0.0       0.0  0.0  0.0;
-            0.866025     0.5       1.0          0.0       0.0       0.0       0.0       0.0  0.0  0.0;
-            0.5          0.866025  0.866025     0.5       1.0       0.0       0.0       0.0  0.0  0.0;
-            2.77556e-16  1.0       0.5          0.866025  0.866025  0.5       1.0       0.0  0.0  0.0;
-           -0.5          0.866025  2.77556e-16  1.0       0.5       0.866025  0.866025  0.5  1.0  0.0
-        ],
-        atol=1e-6
-    ))
+    @test all(
+        isapprox.(
+            X_ϕ1,
+            [
+                1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+                0.866025 0.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+                0.5 0.866025 0.866025 0.5 1.0 0.0 0.0 0.0 0.0 0.0
+                2.77556e-16 1.0 0.5 0.866025 0.866025 0.5 1.0 0.0 0.0 0.0
+                -0.5 0.866025 2.77556e-16 1.0 0.5 0.866025 0.866025 0.5 1.0 0.0
+            ],
+            atol=1e-6,
+        ),
+    )
 end
 
 @testset "Initial State Matrix" begin
-    
-    X1 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, true, Vector{Matrix}(undef, 0))
-    X2 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, true, Vector{Matrix}(undef, 0))
+    X1 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 0, true, true, true, Vector{Matrix}(undef, 0)
+    )
+    X2 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 2, true, true, true, Vector{Matrix}(undef, 0)
+    )
 
     @test X1 == [
         1.0 0.0 1.0 0.0
@@ -278,8 +292,12 @@ end
         1.0 6.0 1.0 0.0
     ]
 
-    X3 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, false, Vector{Matrix}(undef, 0))
-    X4 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, false, Vector{Matrix}(undef, 0))
+    X3 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 0, true, true, false, Vector{Matrix}(undef, 0)
+    )
+    X4 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 2, true, true, false, Vector{Matrix}(undef, 0)
+    )
 
     @test X3 == [
         1.0 0.0
@@ -299,53 +317,69 @@ end
         1.0 6.0
     ]
 
-    X5 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, false, false, Vector{Matrix}(undef, 0))
-    X6 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, false, false, Vector{Matrix}(undef, 0))
+    X5 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 0, true, false, false, Vector{Matrix}(undef, 0)
+    )
+    X6 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 2, true, false, false, Vector{Matrix}(undef, 0)
+    )
 
     @test X5 == ones(5, 1)
     @test X6 == ones(7, 1)
 
-    X7 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, false, true, false, Vector{Matrix}(undef, 0))
-    X8 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, false, true, false, Vector{Matrix}(undef, 0))
+    X7 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 0, false, true, false, Vector{Matrix}(undef, 0)
+    )
+    X8 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 2, false, true, false, Vector{Matrix}(undef, 0)
+    )
 
     @test X7 == [0.0; 1.0; 2.0; 3.0; 4.0][:, :]
     @test X8 == [0.0; 1.0; 2.0; 3.0; 4.0; 5.0; 6.0][:, :]
 
     cycle_matrix = Vector{Matrix}(undef, 1)
     A = 1 * cos(2 * pi / 12)
-    B = 1 * sin(2 * pi / 12) 
+    B = 1 * sin(2 * pi / 12)
     cycle_matrix[1] = [A B; -B A]
     det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0)
     det_matrix2 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 2)
 
-    X9 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, true, true, true, det_matrix1)
-    X10 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, true, true, true, det_matrix2)
+    X9 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 0, true, true, true, det_matrix1
+    )
+    X10 = StateSpaceLearning.create_initial_states_Matrix(
+        5, 2, 2, true, true, true, det_matrix2
+    )
 
-    @test all(isapprox.(
-        X9, 
-        [
-            1.0  0.0  1.0  0.0   1.0          0.0;
-            1.0  1.0  0.0  1.0   0.866025     0.5;
-            1.0  2.0  1.0  0.0   0.5          0.866025;
-            1.0  3.0  0.0  1.0   2.77556e-16  1.0;
-            1.0  4.0  1.0  0.0  -0.5          0.866025;
-        ],
-        atol=1e-6
-    ))
+    @test all(
+        isapprox.(
+            X9,
+            [
+                1.0 0.0 1.0 0.0 1.0 0.0
+                1.0 1.0 0.0 1.0 0.866025 0.5
+                1.0 2.0 1.0 0.0 0.5 0.866025
+                1.0 3.0 0.0 1.0 2.77556e-16 1.0
+                1.0 4.0 1.0 0.0 -0.5 0.866025
+            ],
+            atol=1e-6,
+        ),
+    )
 
-    @test all(isapprox.(
-        X10, 
-        [
-            1.0  0.0  1.0  0.0   1.0          0.0;
-            1.0  1.0  0.0  1.0   0.866025     0.5;
-            1.0  2.0  1.0  0.0   0.5          0.866025;
-            1.0  3.0  0.0  1.0   2.77556e-16  1.0;
-            1.0  4.0  1.0  0.0  -0.5          0.866025;
-            1.0  5.0  0.0  1.0  -0.866025     0.5;
-            1.0  6.0  1.0  0.0  -1.0          5.55112e-16
-        ],
-        atol=1e-6
-    ))
+    @test all(
+        isapprox.(
+            X10,
+            [
+                1.0 0.0 1.0 0.0 1.0 0.0
+                1.0 1.0 0.0 1.0 0.866025 0.5
+                1.0 2.0 1.0 0.0 0.5 0.866025
+                1.0 3.0 0.0 1.0 2.77556e-16 1.0
+                1.0 4.0 1.0 0.0 -0.5 0.866025
+                1.0 5.0 0.0 1.0 -0.866025 0.5
+                1.0 6.0 1.0 0.0 -1.0 5.55112e-16
+            ],
+            atol=1e-6,
+        ),
+    )
 end
 
 @testset "Create X matrix" begin
@@ -422,9 +456,33 @@ end
     counter = 1
     for args in [
         [true, true, true, true, true, true, 2, cycle_matrix, stochastic_cycle, true, 12],
-        [true, true, false, false, false, false, 2, cycle_matrix, stochastic_cycle, true, 12],
+        [
+            true,
+            true,
+            false,
+            false,
+            false,
+            false,
+            2,
+            cycle_matrix,
+            stochastic_cycle,
+            true,
+            12,
+        ],
         [true, true, true, true, false, false, 2, cycle_matrix, stochastic_cycle, true, 12],
-        [true, false, true, true, false, false, 2, cycle_matrix, stochastic_cycle, true, 12],
+        [
+            true,
+            false,
+            true,
+            true,
+            false,
+            false,
+            2,
+            cycle_matrix,
+            stochastic_cycle,
+            true,
+            12,
+        ],
     ]
         args = [x in [0, 1] ? Bool(x) : x for x in args]
         for param in param_combination
@@ -435,9 +493,13 @@ end
                     args..., stochastic_start, Exogenous_X1, param[3], Exogenous_forecast1
                 )
             else
-                X1 = StateSpaceLearning.create_X(args..., stochastic_start, Exogenous_X1, param[3])
+                X1 = StateSpaceLearning.create_X(
+                    args..., stochastic_start, Exogenous_X1, param[3]
+                )
             end
-            X2 = StateSpaceLearning.create_X(args..., stochastic_start, Exogenous_X2, param[3])
+            X2 = StateSpaceLearning.create_X(
+                args..., stochastic_start, Exogenous_X2, param[3]
+            )
             @test size(X1) == size_vec1[counter]
             @test size(X2) == size_vec2[counter]
             counter += 1
@@ -597,10 +659,20 @@ end
         Exogenous_X=Exogenous_X2,
     )
 
-    models = [Basic_Structural, Basic_Structural2, Local_Level, Local_Linear_Trend, Local_Linear_Trend_cycle]
+    models = [
+        Basic_Structural,
+        Basic_Structural2,
+        Local_Level,
+        Local_Linear_Trend,
+        Local_Linear_Trend_cycle,
+    ]
 
     params_vec = [
-        ["ξ", "ζ", "ω_2", "ε"], ["ξ", "ζ", "ω_2", "ω_5", "ε"], ["ξ", "ε"], ["ξ", "ζ", "ε"], ["ξ", "ζ", "ε", "ϕ_1"]
+        ["ξ", "ζ", "ω_2", "ε"],
+        ["ξ", "ζ", "ω_2", "ω_5", "ε"],
+        ["ξ", "ε"],
+        ["ξ", "ζ", "ε"],
+        ["ξ", "ζ", "ε", "ϕ_1"],
     ]
 
     for idx in eachindex(models)
@@ -698,7 +770,11 @@ end
     models = [model1, model2, model3, model4]
 
     keys_vec = [
-        ["ξ", "ζ", "ω_2"], ["ζ", "ω_2"], ["ξ", "ω_2"], ["ξ", "ζ"], ["ξ", "ω_2", "ω_5", "ϕ_1"]
+        ["ξ", "ζ", "ω_2"],
+        ["ζ", "ω_2"],
+        ["ξ", "ω_2"],
+        ["ξ", "ζ"],
+        ["ξ", "ω_2", "ω_5", "ϕ_1"],
     ]
 
     for idx in eachindex(models)
@@ -762,15 +838,17 @@ end
     ]
 
     X4 = StateSpaceLearning.get_innovation_simulation_X(model, innovation4, steps_ahead)
-    @assert all(isapprox.(
-        X4, 
-        [
-            1.0   0.0           0.0   0.0           0.0   0.0;
-            -0.5   0.866025      1.0   0.0           0.0   0.0;
-            -0.5  -0.866025     -0.5   0.866025      1.0   0.0;
-             1.0  -6.10623e-16  -0.5  -0.866025     -0.5   0.866025;
-            -0.5   0.866025      1.0  -6.10623e-16  -0.5  -0.866025;
-        ],
-        atol=1e-6
-    ))
+    @assert all(
+        isapprox.(
+            X4,
+            [
+                1.0 0.0 0.0 0.0 0.0 0.0
+                -0.5 0.866025 1.0 0.0 0.0 0.0
+                -0.5 -0.866025 -0.5 0.866025 1.0 0.0
+                1.0 -6.10623e-16 -0.5 -0.866025 -0.5 0.866025
+                -0.5 0.866025 1.0 -6.10623e-16 -0.5 -0.866025
+            ],
+            atol=1e-6,
+        ),
+    )
 end

From 81c93f21c99b0b8df14ebf0d50e5d2e1346a29fa Mon Sep 17 00:00:00 2001
From: andre_ramos <andreramos@lampsco.com>
Date: Fri, 14 Feb 2025 17:07:51 -0300
Subject: [PATCH 3/4] fix style in paper_tests folder

---
 paper_tests/m4_test/m4_test.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/paper_tests/m4_test/m4_test.jl b/paper_tests/m4_test/m4_test.jl
index cda5c3a..7b18025 100644
--- a/paper_tests/m4_test/m4_test.jl
+++ b/paper_tests/m4_test/m4_test.jl
@@ -140,4 +140,4 @@ end
 
 create_dirs()
 
-main()
+main()

From 133811895df53c3aeffc207fb0273606e49a0276 Mon Sep 17 00:00:00 2001
From: andre_ramos <andreramos@lampsco.com>
Date: Fri, 14 Feb 2025 17:16:13 -0300
Subject: [PATCH 4/4] fix arg int64 bug

---
 src/models/structural_model.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl
index b1d0f5a..0a57f1d 100644
--- a/src/models/structural_model.jl
+++ b/src/models/structural_model.jl
@@ -352,7 +352,7 @@ function create_o_matrix(T::Int, steps_ahead::Int, stochastic_start::Int)::Matri
 end
 
 """
-create_ϕ(X_cycle::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
+create_ϕ(X_cycle::Matrix, T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix
 
     Creates a matrix of innovations ϕ based on the input sizes, and the desired steps ahead (this is necessary for the forecast function).
 
@@ -367,7 +367,7 @@ create_ϕ(X_cycle::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, sto
     - `Matrix`: Matrix of innovations ϕ constructed based on the input sizes.
 """
 function create_ϕ(
-    c_matrix::Matrix, T::Int, steps_ahead::Int64, ζ_ω_threshold::Int, stochastic_start::Int
+    c_matrix::Matrix, T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int
 )::Matrix
     num_cols = 2 * (T - stochastic_start + 1)
     X = Matrix{Float64}(undef, T + steps_ahead, num_cols)