improve tests, fix scaling call (#32)

Author: raphaelsaavedra

docs/src/examples.md

@@ -10,7 +10,7 @@ y = Vector{Float64}(vec(inflow'))
 
 # Specify GAS model: here we use lag 1 for trend characterization and 
 # lag 12 for seasonality characterization
-gas = GAS_Sarima([1, 12], [1, 12], LogNormal, 0.0)
+gas = GAS([1, 12], [1, 12], LogNormal, 0.0)
 
 # Estimate the model via MLE
 estimate!(gas, y)

examples/inflow_lognormal.jl

@@ -48,7 +48,7 @@ inflow = [
 y = Vector{Float64}(vec(inflow'))
 
 # Specify GAS model: here we use lag 1 for trend characterization and lag 12 for seasonality characterization
-gas = GAS_Sarima([1, 12], [1, 12], LogNormal, 0.0)
+gas = GAS([1, 12], [1, 12], LogNormal, 0.0)
 
 # Estimate the model via MLE
 estimate!(gas, y)

src/MLE.jl

@@ -37,7 +37,8 @@ function estimate!(sdm::SDM{D, T}, y::Vector{T};
             push!(psi, optseed.minimizer)
             push!(optseeds, optseed)
             println("seed $i of $nseeds - $(-optseed.minimum)")
-        catch
+        catch err
+            println(err)
             println("seed $i diverged")
         end
     end

src/gas/gas.jl

@@ -5,10 +5,6 @@ mutable struct GAS{D <: Distribution, T <: AbstractFloat} <: SDM{D, T}
     A::Dict{Int, Matrix{T}}
     B::Dict{Int, Matrix{T}}
     scaling::Real
-
-    function GAS{D, T}(ω::Vector{Float64}, A::Dict{Int, Matrix{T}}, B::Dict{Int, Matrix{T}}, scaling::Real) where {D <: Distribution, T <: AbstractFloat}
-        return new(ω, A, B, scaling)
-    end
 end
 
 function create_ω(num_params::Int)

src/score.jl

@@ -17,15 +17,15 @@ function score_tilde(y::T, D::Type{<:Distribution}, param::Vector{T}, param_tild
 end
 
 # Scalings
-function scaling_invsqrt(jac::Matrix{T}, D::Type{Distribution}, param::Vector{T}) where T
+function scaling_invsqrt(jac::AbstractMatrix{T}, D::Type{<:Distribution}, param::Vector{T}) where T
     #TODO improve performace and check if this is right
     # chol = cholesky(fisher_information(dist, param))
     # cholmat = Matrix(chol)
     error("check how to do this properly")
     return jac'*inv(cholmat)*jac
 end
 
-function scaling_inv(jac::Matrix{T}, D::Type{Distribution}, param::Vector{T}) where T
+function scaling_inv(jac::AbstractMatrix{T}, D::Type{<:Distribution}, param::Vector{T}) where T
     #TODO improve performace
     return jac'*inv(fisher_information(D, param))*jac
 end

test/test_distributions.jl

@@ -34,7 +34,7 @@ end
     end
 end
 
-@testset "Log Likelihood" begin
+@testset "Log-likelihood" begin
     for dist in SDM.DISTS
         @testset "$dist" begin
             test_loglik(dist)

test/test_estimate.jl

@@ -1,42 +1,105 @@
-function simulate_GAS_Sarima_1_1(D::Type{Beta}, scaling::Float64)
-    # Create model
-    Random.seed!(13)
-    vec = [0.1 ; 0.1]
-
+function simulate_GAS_1_1(D::Type{<:Distribution}, scaling::Float64, ω::Vector{T}, A::Matrix{T}, 
+                            B::Matrix{T}, seed::Int) where T
+    Random.seed!(seed)
     gas = GAS(1, 1, D, scaling)
 
-    gas.ω = vec
-    gas.A[1] = convert(Matrix{Float64}, Diagonal(5*vec))
-    gas.B[1] = convert(Matrix{Float64}, Diagonal(5*vec))  
+    gas.ω = ω
+    gas.A[1] = A
+    gas.B[1] = B
+    series, param = simulate(gas, 1000)
+
+    return series
+end
+
+function simulate_GAS_1_12(D::Type{<:Distribution}, scaling::Float64)
+    Random.seed!(123)
+    v = [0.1, 0.1]
+    gas = GAS([1, 12], [1, 12], D, scaling)
 
-    # Simulate 1000 observations
-    serie_simulated, param_simulated = simulate(gas, 1000)
+    gas.ω = v
+    gas.A[1]  = convert(Matrix{Float64}, Diagonal(3*v))
+    gas.A[12] = convert(Matrix{Float64}, Diagonal(-3*v))
+    gas.B[1]  = convert(Matrix{Float64}, Diagonal(3*v))
+    gas.B[12] = convert(Matrix{Float64}, Diagonal(-3*v))
 
-    return serie_simulated
+    series, param = simulate(gas, 5000)
+
+    return series
 end
 
-function test_coefficients_GAS_Sarima_1_1(gas::GAS{Beta, T}; atol = 1e-1, rtol = 1e-1) where T
-    @test gas.ω[1] ≈ 0.1 atol = atol rtol = rtol
-    @test gas.ω[2] ≈ 0.1 atol = atol rtol = rtol
-    @test gas.A[1][1, 1] ≈ 0.5 atol = atol rtol = rtol
-    @test gas.A[1][2, 2] ≈ 0.5 atol = atol rtol = rtol
-    @test gas.B[1][1, 1] ≈ 0.5 atol = atol rtol = rtol
-    @test gas.B[1][2, 2] ≈ 0.5 atol = atol rtol = rtol
+function test_coefficients_GAS_1_1(gas::GAS{D, T}, ω::Vector{T}, A::Matrix{T}, B::Matrix{T}; 
+                                    atol = 1e-1, rtol = 1e-1) where {D <: Distribution, T}
+    @test gas.ω[1] ≈ ω[1] atol = atol rtol = rtol
+    @test gas.ω[2] ≈ ω[2] atol = atol rtol = rtol
+    @test gas.A[1][1, 1] ≈ A[1, 1] atol = atol rtol = rtol
+    @test gas.A[1][2, 2] ≈ A[2, 2] atol = atol rtol = rtol
+    @test gas.B[1][1, 1] ≈ B[1, 1] atol = atol rtol = rtol
+    @test gas.B[1][2, 2] ≈ B[2, 2] atol = atol rtol = rtol
     return
 end
 
-function test_estimation_GAS_Sarima_1_1(gas::GAS{D, T}, simulation::Vector{T}) where {D, T}
-    estimate!(gas, simulation; verbose = 1,
-                         opt_method = ScoreDrivenModels.LBFGS(gas, 5))
-
-    test_coefficients_GAS_Sarima_1_1(gas)
-    return 
+function test_coefficients_GAS_1_12(gas::GAS{D, T}; atol = 1e-1, rtol = 1e-1) where {D <: Distribution, T}
+    @test gas.ω[1] ≈ 0.1 atol = atol rtol = rtol
+    @test gas.ω[2] ≈ 0.1 atol = atol rtol = rtol
+    @test gas.A[1][1, 1] ≈ 0.3 atol = atol rtol = rtol
+    @test gas.A[1][2, 2] ≈ 0.3 atol = atol rtol = rtol
+    @test gas.A[12][1, 1] ≈ -0.3 atol = atol rtol = rtol
+    @test gas.A[12][2, 2] ≈ -0.3 atol = atol rtol = rtol
+    @test gas.B[1][1, 1] ≈ 0.3 atol = atol rtol = rtol
+    @test gas.B[1][2, 2] ≈ 0.3 atol = atol rtol = rtol
+    @test gas.B[12][1, 1] ≈ -0.3 atol = atol rtol = rtol
+    @test gas.B[12][2, 2] ≈ -0.3 atol = atol rtol = rtol
+    return
 end
 
 @testset "Estimate" begin
     @testset "Beta" begin
-        simulation = simulate_GAS_Sarima_1_1(Beta, 0.0)
-        gas = GAS(1, 1, Beta, 0.0)
-        test_estimation_GAS_Sarima_1_1(gas, simulation)
+        ω = [0.1, 0.1]
+        A = [0.5 0; 0 0.5]
+        B = [0.5 0; 0 0.5]
+        simulation = simulate_GAS_1_1(Beta, 0.0, ω, A, B, 13)
+        @testset "Estimation by passing number of seeds" begin
+            gas = GAS(1, 1, Beta, 0.0)
+            estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+            test_coefficients_GAS_1_1(gas, ω, A, B)
+        end
+        @testset "Estimation by passing seeds" begin
+            gas = GAS(1, 1, Beta, 0.0)
+            seeds = [[0.1, 0.1, 0.5, 0.5, 0.5, 0.5]]
+            estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, seeds))
+            test_coefficients_GAS_1_1(gas, ω, A, B)
+        end
+    end
+
+    @testset "Lognormal" begin
+        @testset "Scaling = 0.0" begin
+            ω = [0.1, 0.1]
+            A = [0.5 0; 0 0.5]
+            B = [0.5 0; 0 0.5]
+            simulation = simulate_GAS_1_1(LogNormal, 0.0, ω, A, B, 13)
+            gas = GAS(1, 1, LogNormal, 0.0)
+            estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+            test_coefficients_GAS_1_1(gas, ω, A, B)
+        end
+        # @testset "Scaling = 0.5" begin
+        #     gas = GAS(1, 1, LogNormal, 0.5)
+        #     estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+        #     test_coefficients_GAS_1_1(gas)
+        # end
+        # @testset "Scaling = 1.0" begin
+        #     ω = [0.5, 0.5]
+        #     A = [0.1 0; 0 0.1]
+        #     B = [0.1 0; 0 0.1]
+        #     simulation = simulate_GAS_1_1(LogNormal, 1.0, ω, A, B, 123)
+        #     gas = GAS(1, 1, LogNormal, 1.0)
+        #     estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+        #     test_coefficients_GAS_1_1(gas)
+        # end
+        @testset "GAS([1, 12], [1, 12])" begin
+            simulation = simulate_GAS_1_12(LogNormal, 0.0)
+            gas = GAS([1, 12], [1, 12], LogNormal, 0.0)
+            estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+            test_coefficients_GAS_1_12(gas)
+        end
     end
 end