Skip to content

Commit c9f2e4f

Browse files
guilhermebodinguilhermebodin
authored
better logic for the old simulate function (#43)
1 parent 4623d07 commit c9f2e4f

File tree

4 files changed

+102
-68
lines changed

4 files changed

+102
-68
lines changed

‎src/gas/diagnostics.jl‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,12 @@
11
export quantile_residuals, pearson_residuals
22

3-
function quantile_residuals(obs::Vector{T}, gas::GAS{D, T}, initial_params::Matrix{T}) where {D, T}
3+
function quantile_residuals(obs::Vector{T}, gas::GAS{D, T};
4+
initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
45

56
n = length(obs)
67
len_ini_par = length(initial_params)
78
quant_res = Vector{T}(undef, n - len_ini_par)
8-
params_fitted = score_driven_recursion(gas, obs, initial_params)
9+
params_fitted = score_driven_recursion(gas, obs; initial_params = initial_params)
910

1011
for t in axes(params_fitted[len_ini_par + 1:end - 1, :], 1)
1112
dist = update_dist(D, params_fitted[len_ini_par + 1:end - 1, :], t)
@@ -23,12 +24,13 @@ function quantile_residuals(obs::Vector{T}, gas::GAS{D, T}, initial_params::Matr
2324
return quant_res
2425
end
2526

26-
function pearson_residuals(obs::Vector{T}, gas::GAS{D, T}, initial_params::Matrix{T}) where {D, T}
27+
function pearson_residuals(obs::Vector{T}, gas::GAS{D, T};
28+
initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
2729

2830
n = length(obs)
2931
len_ini_par = length(initial_params)
3032
pearson = Vector{T}(undef, n - len_ini_par)
31-
params_fitted = score_driven_recursion(gas, obs, initial_params)
33+
params_fitted = score_driven_recursion(gas, obs; initial_params = initial_params)
3234

3335
for t in axes(params_fitted[len_ini_par + 1:end - 1, :], 1)
3436
dist = update_dist(D, params_fitted[len_ini_par + 1:end - 1, :], t)

‎src/gas/simulate.jl‎

Lines changed: 36 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -1,61 +1,48 @@
1-
export simulate
1+
export simulate, forecast
22

3-
function simulate(gas::GAS{D, T}, n::Int, s::Int) where {D, T}
4-
scenarios_series = Matrix{T}(undef, n, s)
3+
"""
4+
simulate(serie::Vector{T}, gas::GAS{D, T}, N::Int, S::Int, kwargs...) where {D, T}
55
6-
for i in 1:s
7-
serie, params = simulate(gas, n)
8-
scenarios_series[:, i] = serie
6+
Generate scenarios for the future of a time series by updating the GAS recursion `N` times and taking
7+
a sample of the distribution until generate `S` scenarios.
8+
"""
9+
function simulate(serie::Vector{T}, gas::GAS{D, T}, N::Int, S::Int;
10+
initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
11+
# Filter params estimated on the time series
12+
params = score_driven_recursion(gas, serie; initial_params = initial_params)
13+
14+
biggest_lag = number_of_lags(gas)
15+
16+
params_simulation = params[(end - biggest_lag):(end - 1), :]
17+
# Create scenarios matrix
18+
scenarios = Matrix{T}(undef, N, S)
19+
for scenario in 1:s
20+
sim, param = simulate_recursion(gas, N + biggest_lag; initial_params = params_simulation)
21+
scenarios[:, scenario] = sim[biggest_lag + 1:end]
922
end
10-
return scenarios_series
11-
end
1223

13-
function simulate(gas::GAS{D, T}, n::Int) where {D, T}
14-
initial_params = stationary_initial_params(gas)
15-
return simulate(gas, n, initial_params)
24+
return scenarios
1625
end
1726

18-
function simulate(gas::GAS{D, T}, n::Int, initial_params::Matrix{T}) where {D, T}
19-
# Allocations
20-
serie = zeros(n)
21-
n_params = num_params(D)
22-
param = Matrix{T}(undef, n, n_params)
23-
param_tilde = Matrix{T}(undef, n, n_params)
24-
scores_tilde = Matrix{T}(undef, n, n_params)
25-
26-
aux = AuxiliaryLinAlg{T}(n_params)
2727

28-
# Auxiliary Allocation
29-
param_dist = zeros(T, 1, n_params)
28+
"""
29+
forecast(serie::Vector{T}, gas::GAS{D, T}, N::Int; kwargs...) where {D, T}
30+
31+
Forecast future values of a time series by updating the GAS recursion `N` times and
32+
taking the mean of the distribution at each time.
33+
"""
34+
function forecast(serie::Vector{T}, gas::GAS{D, T}, N::Int;
35+
initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
36+
# Filter params estimated on the time series
37+
params = score_driven_recursion(gas, serie; initial_params = initial_params)
3038

3139
biggest_lag = number_of_lags(gas)
3240

33-
# initial_values
34-
for t in 1:biggest_lag
35-
for p in 1:n_params
36-
param[t, p] = initial_params[t, p]
37-
end
38-
link!(param_tilde, D, param, t)
39-
# Sample
40-
updated_dist = update_dist(D, param, t)
41-
serie[t] = sample_observation(updated_dist)
42-
score_tilde!(scores_tilde, serie[t], D, param, aux, gas.scaling, t)
43-
end
44-
45-
update_param_tilde!(param_tilde, gas.ω, gas.A, gas.B, scores_tilde, biggest_lag)
46-
unlink!(param, D, param_tilde, biggest_lag + 1)
47-
updated_dist = update_dist(D, param, biggest_lag + 1)
48-
serie[biggest_lag + 1] = sample_observation(updated_dist)
49-
50-
for i in biggest_lag + 1:n-1
51-
# update step
52-
univariate_score_driven_update!(param, param_tilde, scores_tilde, serie[i], aux, gas, i)
53-
# Sample from the updated distribution
54-
unlink!(param, D, param_tilde, i + 1)
55-
updated_dist = update_dist(D, param, i + 1)
56-
serie[i + 1] = sample_observation(updated_dist)
57-
end
58-
update_param!(param, param_tilde, D, n)
41+
params_simulation = params[(end - biggest_lag):(end - 1), :]
42+
# Create scenarios matrix
43+
forec = Vector{T}(undef, N)
44+
sim, param = simulate_recursion(gas, N + biggest_lag; initial_params = params_simulation, update = mean)
45+
forec = sim[biggest_lag + 1:end]
5946

60-
return serie, param
47+
return forec
6148
end

‎src/gas/univariate_score_driven_recursion.jl‎

Lines changed: 55 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,12 @@
1-
export score_driven_recursion, fitted_mean
1+
export score_driven_recursion, fitted_mean, simulate_recursion
22

33
"""
44
score_driven_recursion(sd_model::SDM, observations::Vector{T}) where T
55
66
start with the stationary params for a
77
"""
8-
function score_driven_recursion(gas::GAS{D, T}, observations::Vector{T}) where {D, T}
9-
initial_params = stationary_initial_params(gas)
10-
return score_driven_recursion(gas, observations, initial_params)
11-
end
12-
13-
function score_driven_recursion(gas::GAS{D, T}, observations::Vector{T}, initial_param::Matrix{T}) where {D, T}
8+
function score_driven_recursion(gas::GAS{D, T}, observations::Vector{T};
9+
initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
1410
@assert gas.scaling in SCALINGS
1511
# Allocations
1612
n = length(observations)
@@ -26,7 +22,7 @@ function score_driven_recursion(gas::GAS{D, T}, observations::Vector{T}, initial
2622
# initial_values
2723
for i in 1:biggest_lag
2824
for p in 1:n_params
29-
param[i, p] = initial_param[i, p]
25+
param[i, p] = initial_params[i, p]
3026
end
3127
link!(param_tilde, D, param, i)
3228
score_tilde!(scores_tilde, observations[i], D, param, aux, gas.scaling, i)
@@ -42,6 +38,54 @@ function score_driven_recursion(gas::GAS{D, T}, observations::Vector{T}, initial
4238
return param
4339
end
4440

41+
42+
function simulate_recursion(gas::GAS{D, T}, n::Int;
43+
initial_params::Matrix{T} = stationary_initial_params(gas),
44+
update::Function = sample_observation) where {D, T}
45+
# Allocations
46+
serie = zeros(n)
47+
n_params = num_params(D)
48+
param = Matrix{T}(undef, n, n_params)
49+
param_tilde = Matrix{T}(undef, n, n_params)
50+
scores_tilde = Matrix{T}(undef, n, n_params)
51+
52+
aux = AuxiliaryLinAlg{T}(n_params)
53+
54+
# Auxiliary Allocation
55+
param_dist = zeros(T, 1, n_params)
56+
57+
biggest_lag = number_of_lags(gas)
58+
59+
# initial_values
60+
for t in 1:biggest_lag
61+
for p in 1:n_params
62+
param[t, p] = initial_params[t, p]
63+
end
64+
link!(param_tilde, D, param, t)
65+
# Sample
66+
updated_dist = update_dist(D, param, t)
67+
serie[t] = update(updated_dist)
68+
score_tilde!(scores_tilde, serie[t], D, param, aux, gas.scaling, t)
69+
end
70+
71+
update_param_tilde!(param_tilde, gas.ω, gas.A, gas.B, scores_tilde, biggest_lag)
72+
unlink!(param, D, param_tilde, biggest_lag + 1)
73+
updated_dist = update_dist(D, param, biggest_lag + 1)
74+
serie[biggest_lag + 1] = update(updated_dist)
75+
76+
for i in biggest_lag + 1:n-1
77+
# update step
78+
univariate_score_driven_update!(param, param_tilde, scores_tilde, serie[i], aux, gas, i)
79+
# Sample from the updated distribution
80+
unlink!(param, D, param_tilde, i + 1)
81+
updated_dist = update_dist(D, param, i + 1)
82+
serie[i + 1] = update(updated_dist)
83+
end
84+
update_param!(param, param_tilde, D, n)
85+
86+
return serie, param
87+
end
88+
4589
function univariate_score_driven_update!(param::Matrix{T}, param_tilde::Matrix{T},
4690
scores_tilde::Matrix{T},
4791
observation::T, aux::AuxiliaryLinAlg{T},
@@ -87,8 +131,9 @@ end
87131
88132
return the fitted mean.. #TODO
89133
"""
90-
function fitted_mean(gas::GAS{D, T}, observations::Vector{T}, initial_params::Matrix{T}) where {D, T}
91-
params_fitted = score_driven_recursion(gas, observations, initial_params)
134+
function fitted_mean(gas::GAS{D, T}, observations::Vector{T};
135+
initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
136+
params_fitted = score_driven_recursion(gas, observations; initial_params = initial_params)
92137
n = size(params_fitted, 1)
93138
fitted_mean = Vector{T}(undef, n)
94139

‎test/utils.jl‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ function simulate_GAS_1_1(D::Type{<:Distribution}, scaling::Float64, ω::Vector{
4949
gas.ω = ω
5050
gas.A[1] = A
5151
gas.B[1] = B
52-
series, param = simulate(gas, 5000)
52+
series, param = simulate_recursion(gas, 5000)
5353

5454
return series
5555
end
@@ -65,7 +65,7 @@ function simulate_GAS_1_12(D::Type{<:Distribution}, scaling::Float64, seed::Int)
6565
gas.B[1] = convert(Matrix{Float64}, Diagonal(3*v))
6666
gas.B[12] = convert(Matrix{Float64}, Diagonal(-3*v))
6767

68-
series, param = simulate(gas, 5000)
68+
series, param = simulate_recursion(gas, 5000)
6969

7070
return series
7171
end
@@ -103,9 +103,9 @@ function normality_quantile_and_pearson_residuals(D, n::Int, lags::Int; seed::In
103103
gas.A[1][[1; 4]] .= 0.2
104104
gas.B[1][[1; 4]] .= 0.2
105105

106-
y, params = simulate(gas, n)
107-
quant_res = quantile_residuals(y, gas, params[1:1, :])
108-
pearson = pearson_residuals(y, gas, params[1:1, :])
106+
y, params = simulate_recursion(gas, n)
107+
quant_res = quantile_residuals(y, gas; initial_params = params[1:1, :])
108+
pearson = pearson_residuals(y, gas, initial_params = params[1:1, :])
109109

110110
# quantile residuals
111111
jb = JarqueBeraTest(quant_res)

0 commit comments

Comments
 (0)