improve examples (#107)

Author: raphaelsaavedra

examples/Project.toml

@@ -8,4 +8,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Plots = "= 1.4.3"
-ScoreDrivenModels = "= 0.1.1"
+ScoreDrivenModels = "= 0.1.1"

examples/ane.jl

@@ -1,17 +1,30 @@
-using ScoreDrivenModels, Plots, DelimitedFiles, Dates, Random
+using Dates, DelimitedFiles, Plots, Random, ScoreDrivenModels
 
+# Define dates and load historical Affluent Natural Energy data
 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);
+y = vec(readdlm("../test/data/ane_northeastern.csv"))
+y_train = y[1:400]
+y_test = y[401:460]
+
+# Set RNG seed to guarantee consistent results
+Random.seed!(123)
+
+# Specify GAS model: a lognormal model with time-varying μ, constant σ, and lags 1 and 12
+gas = Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0; time_varying_params=[1])
+
+# Obtain initial parameters to start the GAS recursion
+initial_params = dynamic_initial_params(y_train, gas)
+
+# Fit specified model to historical data using initial parameters
+f = ScoreDrivenModels.fit!(gas, y_train; initial_params=initial_params)
+
+# Get estimation statistics
 estimation_stats = fit_stats(f)
 
-forec = ScoreDrivenModels.forecast(y_train, gas, 60; S=1_000, initial_params=initial_params)
+# Simulate 1000 future scenarios and obtain the 5% and 95% quantiles in each time period
+forec = forecast(y_train, gas, 60; S=1000, 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)
+# Plot results
+plot(dates[401:460], forec.observation_scenarios, color="grey", width=0.05, label="", ylims=(0, 70))
+plot!(dates[360:460], y[360:460], label="ANE", color="black", xlabel="Months", ylabel="GWmed", legend=:topright)
+plot!(dates[401:460], forec.observation_quantiles, label=["Quantiles" "" ""], color="red", line=:dash)

examples/cpichg.jl

@@ -1,20 +1,39 @@
-using ScoreDrivenModels, DelimitedFiles, Random, Statistics
-Random.seed!(123); 
-y = vec(readdlm("../test/data/cpichg.csv"));
+using DelimitedFiles, Random, ScoreDrivenModels, Statistics
 
-gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2]);
-f = ScoreDrivenModels.fit!(gas, y);
+# Load historical Consumer Price Index data
+y = vec(readdlm("../test/data/cpichg.csv"))
 
-gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2]);
-f = fit!(gas, y; verbose=2);
+# Set RNG seed to guarantee consistent results
+Random.seed!(123)
 
-gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1; 2]);
-f = fit!(gas, y; opt_method=LBFGS(gas, 5));
+# Specify GAS model: a student's t model with location scale transformation
+# (see /src/distributions/non_native_dists.jl in the repository)
+gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2])
 
+# Fit specified model to historical data
+f = fit!(gas, y)
+
+# Next, we show examples using nondefault keyword arguments
+# Note that we need to re-define `gas` since its parameters were populated by `fit!`
+gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2])
+
+# In this example, we set verbose to 2, prompting the package to output optimization results
+f = fit!(gas, y; verbose=2)
+
+# Re-define `gas` once again for another example
+gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params=[1, 2])
+
+# In this example, we set the optimization method to LBFGS with 5 initial points
+f = fit!(gas, y; opt_method=LBFGS(gas, 5))
+
+# Get estimation statistics
 estimation_stats = fit_stats(f)
 
-forec = forecast(y, gas, 12);
+# Perform forecast via simulations for 12 time periods ahead
+forec = forecast(y, gas, 12)
+
+# We can access the parameter forecasts in `forec.parameter_forecasts`
 forec.parameter_forecast
-forec.observation_scenarios
 
-param = score_driven_recursion(gas, y)
+# Similarly, we can access the simulated observation scenarios
+forec.observation_scenarios

examples/garch.jl

@@ -1,6 +1,8 @@
-using ScoreDrivenModels, Statistics, DelimitedFiles
+using DelimitedFiles, ScoreDrivenModels, Statistics
 
-# Overwrite link interface to use IdentityLink for Normal distribution.
+# This is a more advanced example. We want to specify a GARCH model using the GAS framework.
+# To that end, we need to overwrite the default link-unlink interface as it would lead to
+# a different parametrization.
 function ScoreDrivenModels.link!(param_tilde::Matrix{T}, ::Type{Normal}, param::Matrix{T}, t::Int) where T
     param_tilde[t, 1] = link(IdentityLink, param[t, 1])
     param_tilde[t, 2] = link(IdentityLink, param[t, 2])
@@ -17,12 +19,28 @@ function ScoreDrivenModels.jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{Normal
     return
 end
 
-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));
+# Load (TODO: what is this data?) data
+y = vec(readdlm("../test/data/BG96.csv"))
+
+# Set initial parameters as observed mean and variance of the series
+initial_params = [mean(y) var(y)]
+
+# We can set upper and lower bounds for the parameters
+ub = [1.0, 1.0, 0.5, 1.0]
+lb = [-1.0, 0.0, 0.0, 0.5]
+
+# Specify GAS model: a normal model with time-varying σ and lag 1
+gas = Model(1, 1, Normal, 1.0, time_varying_params = [2])
+
+# Set initial point for the optimization
+initial_point = [0.0, 0.5, 0.25, 0.75]
+
+# Fit specified model to historical data
+f = fit!(
+    gas, y;
+    initial_params = initial_params,
+    opt_method = IPNewton(gas, [initial_point]; ub=ub, lb=lb)
+)
+
+# Get estimation statistics
 estimation_stats = fit_stats(f)