GitHub

TVMultiSARMA.jl

Description

This package implements Gibbs sampling for Bayesian inference in time-varying multi-seasonal ARMA (TV-Multi-SARMA) models using the Sequential Monte Carlo (SMC) samplers in SMCsamplers.jl and the Gibbs sampling for dynamic global-local shrinkage priors in DynamicGlobalLocalShrinkage.jl.

Installation

Install from the Julia package manager (via Github) by typing ] in the Julia REPL:

] add git@github.com:mattiasvillani/TVMultiSARMA.jl.git

Multi-seasonal AR model

The TVSAR(p,P) with $p$ regular lags and $P$ seasonal lags with a single seasonality $s$ is

\[\begin{equation*} \phi_t(L)\Phi_t(L^s)y_t = \varepsilon_t , \quad \varepsilon_t \sim N(0,\sigma_t^2) \end{equation*}\]

where $\phi_t(L) = 1 - \phi_{1t} L - \phi_{2t} L^2 - \ldots - \phi_{pt} L^p$ is the regular AR polynomial and $\Phi_t(L^s) = 1 - \Phi_{1t} L^s - \Phi_2 L^{2s} - \ldots - \Phi_{Pt} L^{Ps}$ is the seasonal AR polynomial.

The TVMultiSARMA.jl package allows any number of seasonalities $s_1, s_2, \ldots, s_K$ by including additional seasonal AR polynomials.

The AR parameters $\phi_{jt}$ and $\Phi_{jt}$ can optionally be restricted so that the process is stable at every $t$ by the composite map

\[\boldsymbol{\theta}_t \rightarrow \boldsymbol{r}_t \rightarrow \boldsymbol{\phi}_t, \]

where the unrestricted parameters 𝛉ₜ in ℝᵖ are first mapped to the partial autocorrelations 𝐫ₜ in [−1, 1]ᵖ which are then mapped to the stable AR parameters 𝛟ₜ. The map from 𝛉ₜ to 𝐫ₜ can take many forms, for example, the Monahan transformation

\[ r_{jt} = \frac{\theta_{jt}}{\sqrt{1 + \theta_{jt}^2}}.\]

The unrestricted parameters 𝛉ₜ evolve over time following independent dynamic shrinkage process (DSP) priors. For example, for a single AR parameter $\theta_t$ the DSP prior is

\[\begin{align*} \theta_t &= \theta_{t-1} + \nu_t , \quad \nu_t \sim N(0,\exp(h_t)) \\ h_t &= \mu + \kappa(h_{t-1} -\mu) + \eta_t, \quad \eta_t \sim Z(\alpha,\beta,0,1) \end{align*}\]

where $Z(\alpha,\beta,0,1)$ is the Z-distribution with parameters $\alpha=\beta=1/2$. The DSP prior has a global log-variance $\mu$ that determines the overall degree of time-variation, and a local log-variance component $\eta_t$ that allows for large changes in the parameter innovation variance, changes that can be persistent over time due to the AR(1) structure in $h_t$.

The TVMultiSARMA.jl allows for a stochastic volatility model for the measurement variance $\sigma_t^2$

\[\begin{equation*} \log \sigma_t^2 = \log \sigma_{t-1}^2 + \xi_t, \quad \xi_t \sim N(0,\sigma_\xi^2) \end{equation*}\]

Future versions will include dynamic shrinkage process prior for $\sigma_t^2$.

Limitations and future extensions

The TVMultiSARMA.jl package is limited to time-varying AR models and Bayesian inference using the conditional likelihood. The stochastic volatility model for the measurement errors uses a homoscedastic Gaussian parameter evolation for $\log\sigma_t$. Future versions will extend this by adding:

  • moving average MA and seasonal MA components
  • the exact likelihood
  • global-local shrinkage priors for the stochastic volatility model

The current package is not optimized for speed, and is rather sloppy with memory allocations and type instabilities. Future versions will include speed optimizations.

TVMultiSARMA.ConvertWideMat2Vec — Method
ConvertWideMat2Vec(ϕARpost, p, season)

Converts a T × sum(p) × nIter matrix to a length(season) vector with T × p[l] × nIter elements for the l:th season with lag order p[l]

source
TVMultiSARMA.FindActiveLagsMultiSAR — Method
FindActiveLagsMultiSAR(p, s)

Find the active lags in a general multi-seasonal SAR(p,s) process with seasons in the vector s. The number of lags for each seasonal polynomial is given in the vector p.

source
TVMultiSARMA.GibbsLocalMultiSAR — Function
GibbsLocalMultiSAR(x, modelSettings, priorSettings, algoSettings)

Gibbs sampling from the posterior of the locally stable time-varying multi-seasonal AR(p) model with dynamic shrinkage process prior. Optionally, with stochastic volatility.

source
TVMultiSARMA.MultiSARtoReg — Method
MultiSARtoReg(θ::Vector, p::Vector, s::Vector, activeLags; pacf_map = "monahan", 
    negative_signs = true)

Takes a vector θ of length sum(p) with unrestricted AR/SAR coefficients and returns the non-zero coefficients (determined by activeLags) in the product polynomial Πₗ(1 - ϕₗ₁B^(sₗ) - ϕₗ₂B^(2*sₗ) - ....) = (1 - ϕ̃₁B - ϕ̃₂B² - ....) where the parameters in each AR polynomial is restricted to stability region.

source
TVMultiSARMA.NormalApproxUniformStationary — Method
μ₀, Σ₀ = NormalApproxUniformStationary(p)

Returns mean vector (μ₀) and covariance matrix (Σ₀) in the θ ~ Normal approximation to the uniform prior for the unrestricted θₖ in the multi-seasonal AR with p regular lags and seasonal lags.

For example, p = [2,2,1] is SAR model with two regular lags, two lags at the first seasonal period (e.g. daily cycle) and a single lag on the second seasonal period (e.g. yearly cycle). The approximation minimizes the Hellinger distance to the t and skew-t priors.

Examples

julia> μ₀, Σ₀ = NormalApproxUniformStationary([2, 1])
([0.0, -0.529625, 0.0], [1.0857014809 0.0 0.0; 0.0 0.7363356099999999 0.0; 0.0 0.0 1.0857014809])
source
TVMultiSARMA.PostSpecDensMultiSAR — Method
PostSpecDensMultiSAR(ϕARpost, σₑpost, p, season; ωgrid = 0.01:0.01:π, thinFactor = 100)

Compute posterior draws of the spectral density for a time-varying multi-season SAR model from posterior draws of the AR coefficients (ϕARpost a T×sum(p)×nIter matrix) at every thinFactor iteration.

source
TVMultiSARMA.SARMAasReg — Method
ϕ̂, sₑ = SARMAasReg(x, pFit, season, imposeZeros = true)
Fits a SARMA model as a regression and returns the OLS estimates of the AR coefficients and the standard deviation of the residuals.
source
TVMultiSARMA.SetupARReg — Method
y, Z, T = SetupARReg(x, p)

Sets up the lagged matrix Z and the response vector y for an AR(p) process.

The AR(p) process xt = ϕ₁*x{t-1} + ... + ϕp*x{t-p} + ν_t is transformed into the regression model:

y = Z*ϕ + ν

where y = [x{p+1}, x{p+2}, ..., xT], the t:th row of Z is Z[t,:] = [x{t-1}, x{t-2}, ..., x{t-p}] and ϕ = [ϕ₁, ϕ₂, ..., ϕ_p].

Note: we loose the first p observations in the transformation and the return T is adjusted accordingly to T-p.

Examples

julia> x = simARMA([0.7,-0.2], [0.0], 0.0, 1.0, 100)
julia> y, Z, T = SetupARReg(x, p)
julia> inv(Z'*Z)*Z'y
source
TVMultiSARMA.plotEvolStabilityRegion — Function
plotEvolStabilityRegion(ϕevol, θevol, iter, plotEvolLine = true)

Plot the parameter evoluation in θ-space and ϕ-space (triangle in AR(2) case) up to iteration iter. plotEvolLine = true plots a line connecting the last 10 iterations.

source