1+ @doc raw """
2+ BasicStructuralExplanatory(y::Vector{Fl}, s::Int, X::Matrix{Fl}) where Fl
3+
4+ It is defined by:
5+ ```math
6+ \b egin{gather*}
7+ \b egin{aligned}
8+ y_{t} &= \m u_{t} + \g amma_{t} + \b eta_{t, i}X_{t, i} \v arepsilon_{t} \q uad &\v arepsilon_{t} \s im \m athcal{N}(0, \s igma^2_{\v arepsilon})\\
9+ \m u_{t+1} &= \m u_{t} + \n u_{t} + \x i_{t} \q uad &\x i_{t} \s im \m athcal{N}(0, \s igma^2_{\x i})\\
10+ \n u_{t+1} &= \n u_{t} + \z eta_{t} \q uad &\z eta_{t} \s im \m athcal{N}(0, \s igma^2_{\z eta})\\
11+ \g amma_{t+1} &= -\s um_{j=1}^{s-1} \g amma_{t+1-j} + \o mega_{t} \q uad & \o mega_{t} \s im \m athcal{N}(0, \s igma^2_{\o mega})\\
12+ \b eta_{t+1} &= \b eta_{t}
13+ \e nd{aligned}
14+ \e nd{gather*}
15+ ```
16+
17+ # Example
18+ ```jldoctest
19+ julia> model = BasicStructuralExplanatory(rand(100), 12, rand(100, 2))
20+ BasicStructuralExplanatory model
21+ ```
22+
23+ # References
24+ * Durbin, James, & Siem Jan Koopman. (2012). "Time Series Analysis by State Space Methods: Second Edition." Oxford University Press.
25+ """
26+ mutable struct BasicStructuralExplanatory <: StateSpaceModel
27+ hyperparameters:: HyperParameters
28+ system:: LinearUnivariateTimeVariant
29+ seasonality:: Int
30+ results:: Results
31+ exogenous:: Matrix
32+
33+ function BasicStructuralExplanatory (y:: Vector{Fl} , s:: Int , X:: Matrix{Fl} ) where Fl
34+
35+ @assert length (y) == size (X, 1 )
36+
37+ num_observations = size (X, 1 )
38+ num_exogenous = size (X, 2 )
39+
40+ Z = [vcat ([1 ; 0 ; 1 ; zeros (Fl, s - 2 )], X[t, :]) for t in 1 : num_observations]
41+ T = [[
42+ 1 1 zeros (Fl, 1 , s - 1 ) zeros (Fl, 1 , num_exogenous)
43+ 0 1 zeros (Fl, 1 , s - 1 ) zeros (Fl, 1 , num_exogenous)
44+ 0 0 - ones (Fl, 1 , s - 1 ) zeros (Fl, 1 , num_exogenous)
45+ zeros (Fl, s - 2 , 2 ) Matrix {Fl} (I, s - 2 , s - 2 ) zeros (Fl, s - 2 ) zeros (Fl, 10 , num_exogenous)
46+ zeros (Fl, num_exogenous, 13 ) Matrix {Fl} (I, num_exogenous, num_exogenous)
47+ ] for _ in 1 : num_observations]
48+ R = [[
49+ Matrix {Fl} (I, 3 , 3 )
50+ zeros (Fl, s - 2 , 3 )
51+ zeros (num_exogenous, 3 )
52+ ] for _ in 1 : num_observations]
53+ d = [zero (Fl) for _ in 1 : num_observations]
54+ c = [zeros (Fl, size (T[1 ], 1 )) for _ in 1 : num_observations]
55+ H = [one (Fl) for _ in 1 : num_observations]
56+ Q = [zeros (Fl, 3 , 3 ) for _ in 1 : num_observations]
57+
58+ system = LinearUnivariateTimeVariant {Fl} (y, Z, T, R, d, c, H, Q)
59+
60+ names = [[" sigma2_ε" " sigma2_ξ" " sigma2_ζ" " sigma2_ω" " β_$i " for i in 1 : num_exogenous]]
61+
62+ hyperparameters = HyperParameters {Fl} (names)
63+
64+ return new (hyperparameters, system, s, Results {Fl} (), X)
65+ end
66+ end
67+
68+ function default_filter (model:: BasicStructuralExplanatory )
69+ Fl = typeof_model_elements (model)
70+ steadystate_tol = Fl (1e-5 )
71+ a1 = zeros (Fl, num_states (model))
72+ P1 = Fl (1e6 ) .* Matrix {Fl} (I, num_states (model), num_states (model))
73+ return UnivariateKalmanFilter (a1, P1, num_states (model), steadystate_tol)
74+ end
75+
76+ function initial_hyperparameters! (model:: BasicStructuralExplanatory )
77+ Fl = typeof_model_elements (model)
78+ initial_hyperparameters = Dict {String,Fl} (
79+ " sigma2_ε" => one (Fl),
80+ " sigma2_ξ" => one (Fl),
81+ " sigma2_ζ" => one (Fl),
82+ " sigma2_ω" => one (Fl),
83+ )
84+ initial_exogenous = model. exogenous \ model. system. y
85+ for i in axes (model. exogenous, 2 )
86+ initial_hyperparameters[get_beta_name (model, i)] = initial_exogenous[i]
87+ end
88+ set_initial_hyperparameters! (model, initial_hyperparameters)
89+ return model
90+ end
91+
92+ function get_beta_name (model:: BasicStructuralExplanatory , i:: Int )
93+ return model. hyperparameters. names[i + 4 ]
94+ end
95+
96+ function constrain_hyperparameters! (model:: BasicStructuralExplanatory )
97+ for i in axes (model. exogenous, 2 )
98+ constrain_identity! (model, get_beta_name (model, i))
99+ end
100+ constrain_variance! (model, " sigma2_ε"
101+ constrain_variance! (model, " sigma2_ξ"
102+ constrain_variance! (model, " sigma2_ζ"
103+ constrain_variance! (model, " sigma2_ω"
104+ return model
105+ end
106+
107+ function unconstrain_hyperparameters! (model:: BasicStructuralExplanatory )
108+ for i in axes (model. exogenous, 2 )
109+ unconstrain_identity! (model, get_beta_name (model, i))
110+ end
111+ unconstrain_variance! (model, " sigma2_ε"
112+ unconstrain_variance! (model, " sigma2_ξ"
113+ unconstrain_variance! (model, " sigma2_ζ"
114+ unconstrain_variance! (model, " sigma2_ω"
115+ return model
116+ end
117+
118+ function fill_model_system! (model:: BasicStructuralExplanatory )
119+ H = get_constrained_value (model, " sigma2_ε"
120+ fill_H_in_time (model, H)
121+ for t in 1 : length (model. system. Q)
122+ model. system. Q[t][1 ] = get_constrained_value (model, " sigma2_ξ"
123+ model. system. Q[t][5 ] = get_constrained_value (model, " sigma2_ζ"
124+ model. system. Q[t][end ] = get_constrained_value (model, " sigma2_ω"
125+ end
126+ return model
127+ end
128+
129+ function fill_H_in_time (model:: BasicStructuralExplanatory , H:: Fl ) where Fl
130+ return fill_system_matrice_with_value_in_time (model. system. H, H)
131+ end
132+
133+ function reinstantiate (model:: BasicStructuralExplanatory , y:: Vector{Fl} , X:: Matrix{Fl} ) where Fl
134+ return BasicStructuralExplanatory (y, model. seasonality, X)
135+ end
136+
137+ has_exogenous (:: BasicStructuralExplanatory ) = true
138+
139+ # BasicStructuralExplanatory requires a custom simulation
140+
141+ function simulate_scenarios (
142+ model:: BasicStructuralExplanatory ,
143+ steps_ahead:: Int ,
144+ n_scenarios:: Int ,
145+ new_exogenous:: Matrix{Fl} ;
146+ filter:: KalmanFilter = default_filter (model),
147+ ) where Fl
148+ @assert steps_ahead == size (new_exogenous, 1 ) " new_exogenous must have the same dimension as steps_ahead"
149+ # Query the type of model elements
150+ fo = kalman_filter (model)
151+ last_state = fo. a[end ]
152+ num_series = size (model. system. y, 2 )
153+
154+ scenarios = Array {Fl,3} (undef, steps_ahead, num_series, n_scenarios)
155+ for s in 1 : n_scenarios
156+ scenarios[:, :, s] = simulate (model, last_state, steps_ahead, new_exogenous)
157+ end
158+ return scenarios
159+ end
160+
161+ function simulate_scenarios (
162+ model:: BasicStructuralExplanatory ,
163+ steps_ahead:: Int ,
164+ n_scenarios:: Int ,
165+ new_exogenous:: Array{Fl, 3} ;
166+ filter:: KalmanFilter = default_filter (model),
167+ ) where Fl
168+ @assert steps_ahead == size (new_exogenous, 1 ) " new_exogenous must have the same dimension of steps_ahead"
169+ @assert n_scenarios == size (new_exogenous, 3 ) " new_exogenous must have the same number of scenarios of n_scenarios"
170+ # Query the type of model elements
171+ fo = kalman_filter (model)
172+ last_state = fo. a[end ]
173+ num_series = size (model. system. y, 2 )
174+
175+ scenarios = Array {Fl,3} (undef, steps_ahead, num_series, n_scenarios)
176+ for s in 1 : n_scenarios
177+ scenarios[:, :, s] = simulate (model, last_state, steps_ahead, new_exogenous[:, :, s])
178+ end
179+ return scenarios
180+ end
181+
182+ function simulate (
183+ model:: BasicStructuralExplanatory ,
184+ initial_state:: Vector{Fl} ,
185+ n:: Int ,
186+ new_exogenous:: Matrix{Fl} ;
187+ return_simulated_states:: Bool = false ,
188+ ) where Fl
189+ sys = model. system
190+ m = size (sys. T[1 ], 1 )
191+ y = Vector {Fl} (undef, n)
192+ alpha = Matrix {Fl} (undef, n + 1 , m)
193+ # Sampling errors
194+ chol_H = sqrt (sys. H[1 ])
195+ chol_Q = cholesky (sys. Q[1 ])
196+ standard_ε = randn (n)
197+ standard_η = randn (n + 1 , size (sys. Q[1 ], 1 ))
198+
199+ # The first state of the simulation is the update of a_0
200+ alpha[1 , :] .= initial_state
201+ sys. Z[1 ][14 : end ] .= new_exogenous[1 , :]
202+ y[1 ] = dot (sys. Z[1 ], initial_state) + sys. d[1 ] + chol_H * standard_ε[1 ]
203+ alpha[2 , :] = sys. T[1 ] * initial_state + sys. c[1 ] + sys. R[1 ] * chol_Q. L * standard_η[1 , :]
204+ # Simulate scenarios
205+ for t in 2 : n
206+ sys. Z[t][14 : end ] .= new_exogenous[t, :]
207+ y[t] = dot (sys. Z[t], alpha[t, :]) + sys. d[t] + chol_H * standard_ε[t]
208+ alpha[t + 1 , :] = sys. T[t] * alpha[t, :] + sys. c[t] + sys. R[t] * chol_Q. L * standard_η[t, :]
209+ end
210+
211+ if return_simulated_states
212+ return y, alpha[1 : n, :]
213+ end
214+ return y
215+ end
0 commit comments