Fix idx forecast

Project.toml

@@ -1,7 +1,7 @@
 name = "StateSpaceModels"
 uuid = "99342f36-827c-5390-97c9-d7f9ee765c78"
 authors = ["raphaelsaavedra <raphael.saavedra93@gmail.com>, guilhermebodin <guilherme.b.moraes@gmail.com>, mariohsouto"]
-version = "0.5.6"
+version = "0.5.7"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/forecast.jl

@@ -46,11 +46,11 @@ function forecast(
     for i in 1:steps_ahead
         if isunivariate(model)
             expected_value[i] = [
-                dot(model.system.Z, fo.a[end - steps_ahead + i]) + model.system.d
+                dot(model.system.Z, fo.a[end - steps_ahead + i - 1]) + model.system.d
             ]
         else
             expected_value[i] =
-                model.system.Z * fo.a[end - steps_ahead + i] .+ model.system.d
+                model.system.Z * fo.a[end - steps_ahead + i - 1] .+ model.system.d
         end
         covariance[i] = fo.F[end - steps_ahead + i]
     end
@@ -89,13 +89,13 @@ function forecast(
     for i in 1:steps_ahead
         if isunivariate(model)
             expected_value[i] = [
-                dot(model.system.Z[end - steps_ahead + i], fo.a[end - steps_ahead + i]) +
-                model.system.d[end - steps_ahead + i],
+                dot(model.system.Z[end - steps_ahead + i - 1], fo.a[end - steps_ahead + i - 1]) +
+                model.system.d[end - steps_ahead + i - 1],
             ]
         else
             expected_value[i] =
-                model.system.Z[end - steps_ahead + i] * fo.a[end - steps_ahead + i] +
-                model.system.d[end - steps_ahead + i]
+                model.system.Z[end - steps_ahead + i - 1] * fo.a[end - steps_ahead + i - 1] +
+                model.system.d[end - steps_ahead + i - 1]
         end
         covariance[i] = fo.F[end - steps_ahead + i]
     end

src/systems.jl

@@ -305,9 +305,11 @@ function simulate(
     standard_η = randn(n + 1, size(sys.Q, 1))
 
     # The first state of the simulation is the update of a_0
-    alpha[1, :] = sys.T * initial_state + sys.c + sys.R * chol_Q.L * standard_η[1, :]
+    alpha[1, :] .= initial_state
+    y[1] = dot(sys.Z, initial_state) + sys.d + chol_H * standard_ε[1]
+    alpha[2, :] = sys.T * initial_state + sys.c + sys.R * chol_Q.L * standard_η[1, :]
     # Simulate scenarios
-    for t in 1:n
+    for t in 2:n
         y[t] = dot(sys.Z, alpha[t, :]) + sys.d + chol_H * standard_ε[t]
         alpha[t + 1, :] = sys.T * alpha[t, :] + sys.c + sys.R * chol_Q.L * standard_η[t, :]
     end

test/models/locallevel.jl

@@ -17,6 +17,12 @@
     fit!(model)
     @test loglike(model) ≈ -632.5376 atol = 1e-5 rtol = 1e-5
 
+    forec = forecast(model, 10)
+    @test monotone_forecast_variance(forec)
+    # simulating
+    scenarios = simulate_scenarios(model, 10, 10_000)
+    test_scenarios_adequacy_with_forecast(forec, scenarios)
+
     filter = kalman_filter(model)
 
     # Test that getter functions now work since model has been fitted
@@ -68,6 +74,7 @@
         @test filter.Ptt[t][1] > smoother.V[t][1] # by construction
     end
 
+    # Forecasting of series with missing values
     # forecasting
     forec = forecast(model, 10)
     @test monotone_forecast_variance(forec)

test/models/locallineartrend.jl

@@ -15,6 +15,7 @@
     @test monotone_forecast_variance(forec)
 
     predicted_mean = [
+        5.933257,
         5.897660,
         5.862062,
         5.826465,
@@ -23,8 +24,7 @@
         5.719672,
         5.684074,
         5.648477,
-        5.612879,
-        5.577282,
+        5.612879
     ]
     @test maximum(abs.(predicted_mean .- vcat(forec.expected_value...))) < 1.5e-3
 end

test/models/sarima.jl

@@ -12,6 +12,7 @@
     @test monotone_forecast_variance(forec)
     # Prediction from Pyhton statsmodels
     predicted_mean = [
+        -1.11949821,
         -0.72809028,
         -0.47352952,
         -0.30797034,
@@ -20,16 +21,15 @@
         -0.08472164,
         -0.05510058,
         -0.03583588,
-        -0.02330665,
-        -0.01515799,
+        -0.02330665
     ]
     @test predicted_mean ≈ vcat(forec.expected_value...) atol = 1e-3
     # simualting
-    scenarios = simulate_scenarios(model, 10, 100_000)
+    @time scenarios = simulate_scenarios(model, 10, 30_000)
     # Values are very close to 0.0 so we test with absolute tolerance
     # It attains 1e-3 when we make 10M simulations, which is too much
     # computation for a rather simple test.
-    test_scenarios_adequacy_with_forecast(forec, scenarios; atol=1e-1)
+    test_scenarios_adequacy_with_forecast(forec, scenarios; atol=5e-2)
 
     missing_obs = [6, 16, 26, 36, 46, 56, 66, 72, 73, 74, 75, 76, 86, 96]
     missing_dinternet = copy(dinternet)