Update simulate and select_sample methods

Author: qntwrsm

Project.toml

@@ -23,11 +23,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 
 [compat]
-Distributions = "0.25"
-FillArrays = "0, 1"
-MultivariateStats = "0"
-Statistics = "1"
-StatsAPI = "1"
 julia = "1.9"
 
 [extras]

src/DynamicFactorModels.jl

@@ -17,7 +17,7 @@ using Distributed
 using LinearAlgebra
 using FillArrays
 using Random
-using Distributions
+using Distributions: MvNormal
 
 using StatsAPI: StatisticalModel
 

src/types.jl

@@ -155,6 +155,13 @@ function poly(ε::SpatialMovingAverage)
         return I + Diagonal(spatial(ε)) * weights(ε)
     end
 end
+Base.copy(ε::Simple) = Simple(copy(cov(Σ)))
+function Base.copy(ε::SpatialAutoregression)
+    SpatialAutoregression(copy(spatial(ε)), ε.ρmax, copy(weights(ε)), copy(cov(ε)))
+end
+function Base.copy(ε::SpatialMovingAverage)
+    SpatialMovingAverage(copy(spatial(ε)), ε.ρmax, copy(weights(ε)), copy(cov(ε)))
+end
 
 # factor process
 """
@@ -329,37 +336,37 @@ dynamics(F::AbstractFactorProcess) = F.ϕ
 dynamics(F::UnrestrictedUnitRoot) = I
 dynamics(F::NelsonSiegelUnitRoot) = I
 cov(F::AbstractFactorProcess) = F.Σ
-function copy(F::UnrestrictedStationaryIdentified)
+function Base.copy(F::UnrestrictedStationaryIdentified)
     Λ = copy(loadings(F))
     f = copy(factors(F))
     ϕ = copy(dynamics(F))
 
     return UnrestrictedStationaryIdentified(Λ, f, ϕ)
 end
-function copy(F::UnrestrictedStationaryFull)
+function Base.copy(F::UnrestrictedStationaryFull)
     Λ = copy(loadings(F))
     f = copy(factors(F))
     ϕ = copy(dynamics(F))
     Σ = copy(cov(F))
 
     return UnrestrictedStationaryFull(Λ, f, ϕ, Σ)
 end
-function copy(F::UnrestrictedUnitRoot)
+function Base.copy(F::UnrestrictedUnitRoot)
     Λ = copy(loadings(F))
     f = copy(factors(F))
     Σ = copy(cov(F))
 
     return UnrestrictedUnitRoot(Λ, f, Σ)
 end
-function copy(F::NelsonSiegelStationary)
+function Base.copy(F::NelsonSiegelStationary)
     τ = copy(maturities(F))
     f = copy(factors(F))
     ϕ = copy(dynamics(F))
     Σ = copy(cov(F))
 
     return NelsonSiegelStationary(decay(F), τ, f, ϕ, Σ)
 end
-function copy(F::NelsonSiegelUnitRoot)
+function Base.copy(F::NelsonSiegelUnitRoot)
     τ = copy(maturities(F))
     f = copy(factors(F))
     Σ = copy(cov(F))

src/utilities.jl

@@ -1,8 +1,8 @@
 #=
 utilities.jl
 
-    Provides a collection of utility tools for working with dynamic factor 
-    models, such as simulation, filtering, and smoothing. 
+    Provides a collection of utility tools for working with dynamic factor
+    models, such as simulation, filtering, and smoothing.
 
 @author: Quint Wiersma <q.wiersma@vu.nl>
 
@@ -16,10 +16,10 @@ Select a sample `sample` of the data from the dynamic factor model `model`.
 """
 function select_sample(model::DynamicFactorModel, sample::AbstractUnitRange)
     μ = select_sample(mean(model), sample)
-    ε = select_sample(errors(model), sample)
+    ε = copy(errors(model))
     F = select_sample(process(model), sample)
 
-    return DynamicFactorModel(data(model)[:,sample], μ, ε, F)
+    return DynamicFactorModel(data(model)[:, sample], μ, ε, F)
 end
 
 """
@@ -28,78 +28,82 @@ end
 Select a sample `sample` of the data from the mean model `μ`.
 """
 select_sample(μ::ZeroMean, sample::AbstractUnitRange) = ZeroMean(μ.type, μ.n)
-select_sample(μ::Exogenous, sample::AbstractUnitRange) = Exogenous(regressors(μ)[:,sample], size(slopes(μ), 1))
-
-"""
-    select_sample(ε, sample) -> ε
-
-Select a sample `sample` of the data from the error model `ε`.
-"""
-select_sample(ε::Simple, sample::AbstractUnitRange) = Simple(size(resid(ε), 1), length(sample), type=eltype(resid(ε)))
-select_sample(ε::SpatialAutoregression, sample::AbstractUnitRange) = SpatialAutoregression(size(resid(ε), 1), length(sample), copy(weights(ε)), spatial=length(spatial(ε)) == 1 ? :homo : :hetero, type=eltype(resid(ε)))
-select_sample(ε::SpatialMovingAverage, sample::AbstractUnitRange) = SpatialMovingAverage(size(resid(ε), 1), length(sample), copy(weights(ε)), spatial=length(spatial(ε)) == 1 ? :homo : :hetero, type=eltype(resid(ε)))
+function select_sample(μ::Exogenous, sample::AbstractUnitRange)
+    Exogenous(regressors(μ)[:, sample], size(slopes(μ), 1))
+end
 
 """
     select_sample(F, sample) -> F
 
 Select a sample `sample` of the data from the dynamic factor process `F`.
 """
-select_sample(F::UnrestrictedStationaryIdentified, sample::AbstractUnitRange) = UnrestrictedStationary((size(loadings(F), 1), length(sample), size(F)), dependence=:identified, type=eltype(factors(F)))
-select_sample(F::UnrestrictedStationaryFull, sample::AbstractUnitRange) = UnrestrictedStationary((size(loadings(F), 1), length(sample), size(F)), dependence=:full, type=eltype(factors(F)))
-select_sample(F::UnrestrictedUnitRoot, sample::AbstractUnitRange) = UnrestrictedUnitRoot((size(loadings(F), 1), length(sample), size(F)), type=eltype(factors(F)))
-select_sample(F::NelsonSiegelStationary, sample::AbstractUnitRange) = NelsonSiegelStationary(length(sample), maturities(F), type=eltype(factors(F)))
-select_sample(F::NelsonSiegelUnitRoot, sample::AbstractUnitRange) = NelsonSiegelUnitRoot(length(sample), maturities(F), type=eltype(factors(F))) 
+function select_sample(F::UnrestrictedStationaryIdentified, sample::AbstractUnitRange)
+    UnrestrictedStationary((size(loadings(F), 1), length(sample), size(F)),
+                           dependence = :identified, type = eltype(factors(F)))
+end
+function select_sample(F::UnrestrictedStationaryFull, sample::AbstractUnitRange)
+    UnrestrictedStationary((size(loadings(F), 1), length(sample), size(F)),
+                           dependence = :full, type = eltype(factors(F)))
+end
+function select_sample(F::UnrestrictedUnitRoot, sample::AbstractUnitRange)
+    UnrestrictedUnitRoot((size(loadings(F), 1), length(sample), size(F)),
+                         type = eltype(factors(F)))
+end
+function select_sample(F::NelsonSiegelStationary, sample::AbstractUnitRange)
+    NelsonSiegelStationary(length(sample), maturities(F), type = eltype(factors(F)))
+end
+function select_sample(F::NelsonSiegelUnitRoot, sample::AbstractUnitRange)
+    NelsonSiegelUnitRoot(length(sample), maturities(F), type = eltype(factors(F)))
+end
 
 """
-    simulate(F; burn=100, rng=Xoshiro()) -> sim
+    simulate(F, S; rng = Xoshiro()) -> f
 
-Simulate the dynamic factors from the dynamic factor process `F`
-with a burn-in period of `burn` using the random number generator `rng`.
+Simulate the dynamic factors from the dynamic factor process `F` `S` times using the random
+number generator `rng`.
 """
-function simulate(F::AbstractFactorProcess; burn::Integer=100, rng::AbstractRNG=Xoshiro())
-    # burn-in
-    f_prev = rand(rng, dist(F))
-    f_next = similar(f_prev)
-    for _ = 1:burn
-        f_next .= dynamics(F) * f_prev + rand(rng, dist(F))
-        f_prev .= f_next
-    end
-
-    # simulate data
-    F_sim = copy(F)
-    for (t, ft) ∈ pairs(eachcol(factors(F_sim)))
-        if t == 1
-            ft .= f_prev
+function simulate(F::AbstractFactorProcess, S::Integer; rng::AbstractRNG = Xoshiro())
+    R = size(loadings(F), 2)
+    f = similar(factors(F), R, S)
+    dist = MvNormal(cov(F))
+    for (s, fs) in pairs(eachcol(f))
+        if s == 1
+            # initial condition
+            fs .= rand(rng, dist)
         else
-            ft .= dynamics(F) * factors(F_sim)[:,t-1] + rand(rng, dist(F))
+            fs .= dynamics(F) * f[:, s - 1] + rand(rng, dist)
         end
+
     end
 
-    return F_sim
+    return f
 end
 
 """
-    simulate(ε; rng=Xoshiro()) -> sim
+    simulate(ε, S; rng = Xoshiro()) -> e
 
-Simulate from the error distribution `ε` using the random number generator 
-`rng`.
+Simulate from the error distribution `ε` `S` times using the random number generator `rng`.
 """
-simulate(ε::Simple; rng::AbstractRNG=Xoshiro()) = Simple(rand(rng, dist(ε), size(resid(ε), 2)), MvNormal(Diagonal(var(ε))))
-function simulate(ε::SpatialAutoregression; rng::AbstractRNG=Xoshiro())
-    e_sim = similar(resid(ε))
-    for et ∈ eachcol(e_sim)
-        et .= poly(ε) \ rand(rng, dist(ε))
+function simulate(ε::Simple, S::Integer; rng::AbstractRNG = Xoshiro())
+    rand(rng, MvNormal(cov(ε)), S)
+end
+function simulate(ε::SpatialAutoregression, S::Integer; rng::AbstractRNG = Xoshiro())
+    e = rand(rng, MvNormal(cov(ε)), S)
+    G = poly(ε)
+    for es in eachcol(e)
+        es .= G \ es
     end
 
-    return SpatialAutoregression(e_sim, MvNormal(Diagonal(var(ε))), copy(spatial(ε)), ε.ρ_max, copy(weights(ε)))
+    return e
 end
-function simulate(ε::SpatialMovingAverage; rng::AbstractRNG=Xoshiro())
-    e_sim = similar(resid(ε))
-    for et ∈ eachcol(e_sim)
-        mul!(et, poly(ε), rand(rng, dist(ε)))
+function simulate(ε::SpatialMovingAverage, S::Integer; rng::AbstractRNG = Xoshiro())
+    e = rand(rng, MvNormal(cov(ε)), S)
+    G = poly(ε)
+    for es in eachcol(e)
+        es .= G * es
     end
 
-    return SpatialMovingAverage(e_sim, MvNormal(Diagonal(var(ε))), copy(spatial(ε)), ε.ρ_max, copy(weights(ε)))
+    return e
 end
 
 """
@@ -110,15 +114,15 @@ following the approach of Jungbacker and Koopman (2015).
 """
 function collapse(model::DynamicFactorModel)
     # active factors
-    active = [!all(iszero, λ) for λ ∈ eachcol(loadings(model))]
+    active = [!all(iszero, λ) for λ in eachcol(loadings(model))]
 
     # collapsing matrix
     H = cov(model)
     if all(active)
         C = (H \ loadings(model))' * loadings(model)
         Z_basis = (C \ loadings(model)')'
     else
-        Z_basis = loadings(model)[:,active]
+        Z_basis = loadings(model)[:, active]
     end
     A_low = (H \ Z_basis)'
 
@@ -136,7 +140,7 @@ function state_space(model::DynamicFactorModel)
     Ty = eltype(data(model))
 
     # active factors
-    active = [!all(iszero, λ) for λ ∈ eachcol(loadings(model))]
+    active = [!all(iszero, λ) for λ in eachcol(loadings(model))]
 
     # collapsing matrix
     if any(active)
@@ -182,7 +186,7 @@ function _filter(y, Z, d, H, T, Q, a1, P1)
     P[1] = P1
 
     # filter
-    for (t, yt) ∈ pairs(eachcol(y))
+    for (t, yt) in pairs(eachcol(y))
         # forecast error
         v[t] = yt - Z * a[t] - view(d, :, t)
         F[t] = Z * P[t] * Z' + H
@@ -192,8 +196,8 @@ function _filter(y, Z, d, H, T, Q, a1, P1)
 
         # prediction
         if t < n
-            a[t+1] = T * a[t] + K[t] * v[t]
-            P[t+1] = T * P[t] * (T - K[t] * Z)' + Q
+            a[t + 1] = T * a[t] + K[t] * v[t]
+            P[t + 1] = T * P[t] * (T - K[t] * Z)' + Q
         end
     end
 
@@ -242,7 +246,7 @@ function smoother(model::DynamicFactorModel)
     L = similar(P1)
 
     # smoother
-    for t ∈ reverse(eachindex(a))
+    for t in reverse(eachindex(a))
         L .= T - K[t] * Z
 
         # backward recursion
@@ -252,7 +256,7 @@ function smoother(model::DynamicFactorModel)
         # smoothing
         α̂[t] = a[t] + P[t] * r
         V[t] = P[t] - P[t] * N * P[t]
-        t > 1 && (Γ[t-1] = I - P[t] * N)
+        t > 1 && (Γ[t - 1] = I - P[t] * N)
         t < length(a) && (Γ[t] *= L * P[t])
     end
 
@@ -265,4 +269,4 @@ end
 Forecast the mean model `μ` `periods` ahead.
 """
 forecast(μ::AbstractMeanSpecification, periods::Integer) = error("forecast for $(Base.typename(typeof(μ)).wrapper) not implemented.")
-forecast(μ::ZeroMean, periods::Integer) = Zeros(μ.type, μ.n, periods)
+forecast(μ::ZeroMean, periods::Integer) = Zeros(μ.type, μ.n, periods)