New examples (#106)

* update examples

Co-authored-by: Raphael Saavedra <[email protected]>

Author: guilhermebodin

examples/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ScoreDrivenModels = "4a87933e-d659-11e9-0e65-7f40dedd4a3a"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[compat]
+Plots = "= 1.4.3"
+ScoreDrivenModels = "= 0.1.1"

examples/ane.jl

@@ -1,26 +1,17 @@
-using ScoreDrivenModels, Plots
-
-# Convert data to vector
-y = Vector{Float64}(vec(readdlm("./test/data/ena_northeastern.csv")'))
-y_train = y[1:400]
-y_test = y[401:end]
-
-# Specify model: here we use lag 1 for trend characterization and lag 12 for seasonality characterization
-gas = Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0; time_varying_params = [1])
-
-# Define initial_params with
-initial_params = dynamic_initial_params(y_train, gas)
-
-# Estimate the model via MLE
-f = fit!(gas, y_train; initial_params = initial_params)
-
-# Compare observations and in-sample estimates
-plot(y_train, label = "In-sample ANE")
-
-# Forecasts with 95% confidence interval
-forecast = forecast_quantiles(y_train, gas, 60; S = 1_000, initial_params = initial_params)
-
-plot(forecast.scenarios, color = "grey", width = 0.05, label = "")
-plot!(y[360:460], label = "ANE", color = "black", xlabel = "Months", ylabel = "GWmed", legend = :topright)
-plot!(forecast.quantiles, label = ["Quantiles" "" ""], color = "red", line = :dash)
-
+using ScoreDrivenModels, Plots, DelimitedFiles, Dates, Random
+
+dates = collect(Date(1961):Month(1):Date(2000, 12))
+Random.seed!(123);
+y = vec(readdlm("../test/data/ane_northeastern.csv"));
+y_train = y[1:400];
+gas = Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0; time_varying_params=[1]);
+initial_params = dynamic_initial_params(y_train, gas);
+f = ScoreDrivenModels.fit!(gas, y_train; initial_params=initial_params);
+estimation_stats = fit_stats(f)
+
+forec = ScoreDrivenModels.forecast(y_train, gas, 60; S=1_000, initial_params=initial_params)
+
+y_test = y[401:460]
+p2 = plot(dates[401:460], forec.observation_scenarios, color="grey", width=0.05, label="", ylims=(0, 70))
+plot!(p2, dates[360:460], y[360:460], label="ANE", color="black", xlabel="Months", ylabel="GWmed", legend=:topright)
+plot!(p2, dates[401:460], forec.observation_quantiles, label=["Quantiles" "" ""], color="red", line=:dash)

examples/cpichg.jl

@@ -1,11 +1,20 @@
-using ScoreDrivenModels, DelimitedFiles
+using ScoreDrivenModels, DelimitedFiles, Random, Statistics
+Random.seed!(123); 
+y = vec(readdlm("../test/data/cpichg.csv"));
 
-y = readdlm("./test/data/cpichg.csv")[:]
+gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2]);
+f = ScoreDrivenModels.fit!(gas, y);
 
-# Define a model where only \mu and \sigma^2 is varying over time
-gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params = [1; 2])
+gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2]);
+f = fit!(gas, y; verbose=2);
 
-# Estimate the model
-f = fit!(gas, y)
+gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1; 2]);
+f = fit!(gas, y; opt_method=LBFGS(gas, 5));
 
-estimation_stats = fit_stats(f)
+estimation_stats = fit_stats(f)
+
+forec = forecast(y, gas, 12);
+forec.parameter_forecast
+forec.observation_scenarios
+
+param = score_driven_recursion(gas, y)

examples/garch.jl

@@ -17,42 +17,12 @@ function ScoreDrivenModels.jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{Normal
     return
 end
 
-# Data from [Bollerslev and Ghysels (JBES 1996)](https://doi.org/10.2307/1392425). Obtained from [ARCHModels.jl](https://github.com/s-broda/ARCHModels.jl)
-y = readdlm("./test/data/BG96.csv")[:]
-
-# Evaluate some approximate initial params
-initial_params = [mean(y) var(y)]
-
-# There are 2 possibilities to estimate the model.
-# The first one is to understand that in the GARCH estimation process usually one
-# defines that the parameters ω, α_0, α_1 and β_1 are constrained. The problem is
-# that the equivalent GAS model relies on
-# ω_1 = ω
-# ω_2 = α_0
-# A_1 = α_1
-# B_1 - A_1 = β_1
-# Here there is no easy way of adding a constraint B_1 - A_1 >= 0. An alternative is to define
-# bounds assuring the upper bound of A_1 is smaller than the lower bound for B_1.
-
-# Define upper bounds and lower bounds for parameters ω_1, ω_2, A_1, B_1
-ub = [1.0; 1.0; 0.5; 1.0]
-lb = [-1.0; 0.0; 0.0; 0.5]
-
-# Define a model where only \sigma^2 is varying over time
-gas = Model(1, 1, Normal, 1.0, time_varying_params = [2])
-
-# Give an initial_point in the interior of the bounds.
-initial_point = [0.0; 0.5; 0.25; 0.75]
-
-# Estimate the model
-f_ip_newton = fit(gas, y; initial_params = initial_params,
-                  opt_method = IPNewton(gas, [initial_point]; ub = ub, lb = lb))
-
-estimation_stats_ip_newton = fit_stats(f_ip_newton)
-
-# Another way to estimate the model is to rely on luck and give many random initial_points
-# for a non constrained optimization
-f_nelder_mead = fit(gas, y; initial_params = initial_params,
-                  opt_method = NelderMead(gas, 100))
-
-estimation_stats_nelder_mead = fit_stats(f_nelder_mead)
+y = vec(readdlm("../test/data/BG96.csv"));
+initial_params = [mean(y) var(y)];
+ub = [1.0, 1.0, 0.5, 1.0];
+lb = [-1.0, 0.0, 0.0, 0.5];
+gas = Model(1, 1, Normal, 1.0, time_varying_params = [2]);
+initial_point = [0.0, 0.5, 0.25, 0.75];
+f = fit!(gas, y; initial_params = initial_params,
+                opt_method = IPNewton(gas, [initial_point]; ub = ub, lb = lb));
+estimation_stats = fit_stats(f)