Package {wARMASVp}

Type:	Package
Title:	Winsorized ARMA Estimation for Higher-Order Stochastic Volatility Models
Version:	0.2.0
Description:	Estimation, simulation, hypothesis testing, AR-order selection, and forecasting for univariate higher-order stochastic volatility SV(p) models. Supports Gaussian, Student-t, and Generalized Error Distribution (GED) innovations, with optional leverage effects. Estimation uses closed-form Winsorized ARMA-SV (W-ARMA-SV) moment-based methods that avoid numerical optimization. Hypothesis testing includes Local Monte Carlo (LMC) and Maximized Monte Carlo (MMC) procedures for leverage effects, heavy tails, and autoregressive order. AR-order selection is also available via information criteria (BIC/AIC) using the Kalman-filter quasi-likelihood and the Hannan-Rissanen ARMA residual variance. Forecasting is based on Kalman filtering and smoothing. See Ahsan and Dufour (2021) <doi:10.1016/j.jeconom.2021.03.008>, Ahsan, Dufour, and Rodriguez-Rondon (2025) <doi:10.1111/jtsa.12851>, and Ahsan, Dufour, and Rodriguez-Rondon (2026) <doi:10.34989/swp-2026-8> for details.
License:	GPL (≥ 3)
URL:	https://github.com/roga11/wARMASVp
BugReports:	https://github.com/roga11/wARMASVp/issues
Encoding:	UTF-8
Imports:	Rcpp (≥ 1.0.0), gsignal, stats
Suggests:	pso, GenSA, testthat (≥ 3.0.0), knitr, rmarkdown
LinkingTo:	Rcpp, RcppArmadillo
RoxygenNote:	7.3.3
VignetteBuilder:	knitr
NeedsCompilation:	yes
Packaged:	2026-05-15 12:58:15 UTC; gabrielrodriguez
Author:	Gabriel Rodriguez-Rondon [aut, cre], Md. Nazmul Ahsan [aut], Jean-Marie Dufour [aut]
Maintainer:	Gabriel Rodriguez-Rondon <gabriel.rodriguezrondon@mail.mcgill.ca>
Repository:	CRAN
Date/Publication:	2026-05-15 15:50:01 UTC

wARMASVp: Winsorized ARMA Estimation for Higher-Order Stochastic Volatility Models

Description

Estimation, simulation, hypothesis testing, and forecasting for univariate higher-order stochastic volatility SV(p) models. Supports Gaussian, Student-t, and GED innovations with optional leverage effects.

Details

The main user-facing functions are:

• svp – Estimate SV(p) model with optional leverage

• svpSE – Simulation-based standard errors

• sim_svp – Simulate SV(p) processes

• filter_svp – Kalman/mixture/particle filtering

• forecast_svp – Multi-step volatility forecasts

• lmc_ar, mmc_ar – AR order tests

• lmc_lev, mmc_lev – Leverage tests

• lmc_t, mmc_t – Student-t tail tests

• lmc_ged, mmc_ged – GED tail tests

Author(s)

Maintainer: Gabriel Rodriguez-Rondon gabriel.rodriguezrondon@mail.mcgill.ca (ORCID)

Authors:

• Md. Nazmul Ahsan

• Jean-Marie Dufour

References

Ahsan, M. N. and Dufour, J.-M. (2021). Simple estimators and inference for higher-order stochastic volatility models. Journal of Econometrics, 224(1), 181-197. doi:10.1016/j.jeconom.2021.03.008

Ahsan, M. N., Dufour, J.-M., and Rodriguez-Rondon, G. (2025). Estimation and inference for higher-order stochastic volatility models with leverage. Journal of Time Series Analysis, 46(6), 1064-1084. doi:10.1111/jtsa.12851

Ahsan, M. N., Dufour, J.-M., and Rodriguez-Rondon, G. (2026). Estimation and inference for stochastic volatility models with heavy-tailed distributions. Bank of Canada Staff Working Paper 2026-8. doi:10.34989/swp-2026-8

See Also

Useful links:

• https://github.com/roga11/wARMASVp

• Report bugs at https://github.com/roga11/wARMASVp/issues

Add leverage estimation to a fitted SV(p) model (post-processing) Works for all error distributions: Gaussian, Student-t, GED. - Gaussian: closed form (C_F = 1) - Student-t: closed form with C_t(nu) correction - GED: exact 1D root-finding via uniroot + Gauss-Hermite quadrature

Description

Add leverage estimation to a fitted SV(p) model (post-processing) Works for all error distributions: Gaussian, Student-t, GED. - Gaussian: closed form (C_F = 1) - Student-t: closed form with C_t(nu) correction - GED: exact 1D root-finding via uniroot + Gauss-Hermite quadrature

Usage

.add_leverage(out, y, p, rho_type, del, trunc_lev, wDecay, errorType)

Arguments

Value

Updated model object with leverage fields added

Compute EH cross-moment for leverage estimation

Description

EH = demeaned sample cross-moment (1/(T-2)) \sum(|y_t| - \bar{|y|})(y_{t-1} - \bar{y}).

Usage

.compute_EH(y, rho_type = "pearson")

Arguments

Value

EH value

Gauss-Hermite quadrature nodes and weights for N(0,1) integration Computes nodes z_i and weights w_i such that sum(w_i * f(z_i)) approximates E[f(Z)] where Z ~ N(0,1). Uses the Golub-Welsch algorithm.

Description

Gauss-Hermite quadrature nodes and weights for N(0,1) integration Computes nodes z_i and weights w_i such that sum(w_i * f(z_i)) approximates E[f(Z)] where Z ~ N(0,1). Uses the Golub-Welsch algorithm.

Usage

.gauss_hermite_normal(n = 200L)

Arguments

Value

List with components nodes and weights

Get cached Gauss-Hermite nodes/weights for N(0,1)

Description

Get cached Gauss-Hermite nodes/weights for N(0,1)

Usage

.get_gh(n = 200L)

Arguments

Value

List with nodes and weights

GMM moments for SVL(p)-GED with fixed A matrix (p+4 conditions) Uses exact GED leverage moment.

Description

GMM moments for SVL(p)-GED with fixed A matrix (p+4 conditions) Uses exact GED leverage moment.

Usage

LRT_moment_lev_ged(y, mdl_out, Amat, rho_type, del = 1e-10)

GMM moments for SVL(p)-GED with HAC weighting (p+4 conditions)

Description

GMM moments for SVL(p)-GED with HAC weighting (p+4 conditions)

Usage

LRT_moment_lev_ged_Amat(y, mdl_out, rho_type, del = 1e-10, Bartlett = TRUE)

GMM moments for SVL(p)-Student-t with fixed A matrix (p+4 conditions)

Description

GMM moments for SVL(p)-Student-t with fixed A matrix (p+4 conditions)

Usage

LRT_moment_lev_t(y, mdl_out, Amat, rho_type, del = 1e-10)

GMM moments for SVL(p)-Student-t with HAC weighting (p+4 conditions)

Description

GMM moments for SVL(p)-Student-t with HAC weighting (p+4 conditions)

Usage

LRT_moment_lev_t_Amat(y, mdl_out, rho_type, del = 1e-10, Bartlett = TRUE)

Construct Companion Matrix for AR(p) Process

Description

Builds the companion matrix representation for an AR(p) process, which is used in the state-space formulation of SV(p) models.

Usage

companionMat(phi, p, q)

Arguments

Value

A (q*p) x (q*p) companion matrix.

Approximate correction factor C_g(nu) for GED leverage (diagnostic use only)

Description

C_g(\nu) = E[|u|] \cdot \textrm{Cov}(z, F_{GED}^{-1}(\Phi(z))) / \sqrt{2/\pi}. This is a first-order approximation; estimation uses the exact implicit equation.

Usage

correction_factor_ged_approx(nu, n_nodes = 200L)

Arguments

Value

Approximate C_g(nu)

Correction factor C_t(nu) for Student-t leverage estimation

Description

Under scale mixture u_t = z_t \lambda_t^{-1/2}, C_t(\nu) = [E[\lambda^{-1/2}]]^2 = (\nu/2) [\Gamma((\nu-1)/2) / \Gamma(\nu/2)]^2. Exact, parameter-free.

Usage

correction_factor_t(nu)

Arguments

Value

C_t(nu), always > 1 for finite nu (approaches 1 as nu -> Inf)

Density of Centered log-GED^2 Measurement Noise

Description

Density of Centered log-GED^2 Measurement Noise

Usage

density_eps_ged(y, nu)

Arguments

Value

Density values.

Density of Centered log-F(1,nu) Measurement Noise (Student-t)

Description

Density of Centered log-F(1,nu) Measurement Noise (Student-t)

Usage

density_eps_t(y, nu)

Arguments

Value

Density values.

Filter Latent Volatility from an Estimated SV(p) Model

Description

Applies Kalman filtering (corrected or Gaussian mixture) and RTS smoothing to extract the latent log-volatility process from an estimated SV(p) model.

Usage

filter_svp(
  object,
  method = c("corrected", "mixture", "particle"),
  proxy = c("bayes_optimal", "u"),
  K = 7,
  M = 1000,
  seed = 42,
  del = 1e-10
)

Arguments

Value

An object of class "svp_filter", a list containing:

w_filtered

Filtered log-volatility (T-vector).

w_smoothed

Smoothed log-volatility (T-vector).

zt

Filtered standardized residuals.

zt_smoothed

Smoothed standardized residuals.

P_filtered

Filtered MSE of first state component.

P_predicted

Predicted MSE of first state component.

xi_filtered

Full filtered state vectors (p x T matrix).

xi_smoothed

Full smoothed state vectors (p x T matrix).

loglik

Approximate log-likelihood.

method

Filter method used.

model

The input model object.

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2)$y
fit <- svp(y, p = 1)
filt <- filter_svp(fit)
plot(filt$w_smoothed, type = "l")

Fit K-Component Gaussian Mixture to Measurement Noise Density

Description

Uses EM algorithm to approximate the measurement noise density with a Gaussian mixture. For Gaussian SV, returns the pre-computed KSC (1998) 7-component table.

Usage

fit_ksc_mixture(
  distribution = c("gaussian", "student_t", "ged"),
  nu = NULL,
  K = 7,
  n_sample = 10000,
  max_iter = 500,
  tol = 1e-08,
  seed = 42
)

Arguments

Value

List with weights, means, vars, KL_div.

Multi-Step Ahead Volatility Forecast

Description

Applies Kalman filtering/smoothing to an estimated SV(p) model and produces multi-step ahead volatility forecasts with uncertainty quantification.

Usage

forecast_svp(
  object,
  H = 1,
  output = c("log-variance", "variance", "volatility"),
  filter_method = "corrected",
  proxy = c("bayes_optimal", "u"),
  K = 7,
  M = 1000,
  seed = 42,
  del = 1e-10
)

Arguments

Value

An object of class "svp_forecast", a list containing:

w_forecasted

Primary forecast (scale determined by output).

log_var_forecast

Log-volatility forecasts w_{T+h|T}.

var_forecast

Conditional variance forecasts \sigma^2_{T+h|T}.

vol_forecast

Conditional volatility forecasts \sigma_{T+h|T}.

P_forecast

Forecast MSE P_{T+h|T} for each horizon.

w_estimated

Filtered log-volatility.

w_smoothed

Smoothed log-volatility.

zt

Filtered standardized residuals.

zt_smoothed

Smoothed standardized residuals.

ys

Demeaned log-squared returns.

mdl

The estimated model object.

H

The forecast horizon.

output

The chosen output scale.

filter_output

The "svp_filter" object from filtering.

Examples

sim <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2,
               leverage = TRUE, rho = -0.3)
fit <- svp(sim$y, p = 1, leverage = TRUE)
fc <- forecast_svp(fit, H = 10)
plot(fc)

E[|u|] for standardized GED(nu) with Var = 1 Closed form: sqrt(Gamma(1/nu)/Gamma(3/nu)) * Gamma(2/nu) / Gamma(1/nu)

Description

E[|u|] for standardized GED(nu) with Var = 1 Closed form: sqrt(Gamma(1/nu)/Gamma(3/nu)) * Gamma(2/nu) / Gamma(1/nu)

Usage

ged_E_abs_u(nu)

Arguments

Value

E[|u|]

LMC Test for AR Order in SV(p) Models

Description

Performs a Local Monte Carlo (LMC) test of the null hypothesis H_0: \phi_{p_0+1} = \cdots = \phi_p = 0 (i.e., that an SV(p_0) model is sufficient against an SV(p) alternative).

Usage

lmc_ar(
  y,
  p_null,
  p_alt,
  J = 10,
  N = 99,
  burnin = 500,
  del = 1e-10,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  errorType = "Gaussian",
  sigvMethod = "factored",
  logNu = TRUE,
  winsorize_eps = 0
)

Arguments

Details

When Bartlett = FALSE (default), the test statistic is T \sum_{j=p_0+1}^{p} \hat\phi_j^2, i.e., the sum of squared extra AR coefficients scaled by sample size.

When Bartlett = TRUE, the test statistic is based on the GMM-LRT approach with a Bartlett kernel HAC weighting matrix: S = T \times (M_{H_0} - M_{H_1}), where M denotes the GMM criterion evaluated at the null and alternative estimates. Both the observed and simulated test statistics are capped at 1e-10 when negative; a negative observed statistic raises a warning (it indicates strong evidence in favour of the null, since the alternative does not improve the GMM criterion).

Value

An object of class "svp_test", a list containing:

s0

Test statistic from observed data (capped at 1e-10 if negative).

sN

Simulated null distribution (vector of length N).

pval

Monte Carlo p-value.

test_type

Character string identifying the test.

null_param

Name of the parameter(s) tested.

null_value

Value(s) under the null hypothesis.

errorType

Error distribution used.

call

The matched call.

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2)$y
test <- lmc_ar(y, p_null = 1, p_alt = 2, J = 10, N = 49)
print(test)

LMC Test for GED Shape Parameter

Description

Performs a Local Monte Carlo (LMC) test of the null hypothesis H_0: \nu = \nu_0 for the shape parameter in an SV(p) model with GED errors. Testing \nu_0 = 2 corresponds to testing normality.

Usage

lmc_ged(
  y,
  p = 1,
  J = 10,
  N = 99,
  nu_null,
  burnin = 500,
  del = 1e-10,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  direction = c("two-sided", "less", "greater"),
  sigvMethod = "factored",
  winsorize_eps = 0
)

Arguments

Value

An object of class "svp_test".

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2, errorType = "GED", nu = 1.5)$y
test <- lmc_ged(y, p = 1, J = 10, N = 49, nu_null = 2)
print(test)

LMC Test for Leverage in SV(p) Models

Description

Performs a Local Monte Carlo (LMC) test of the null hypothesis H_0: \rho = \rho_0 (typically \rho_0 = 0, i.e., no leverage) using a GMM likelihood-ratio type statistic.

Usage

lmc_lev(
  y,
  p = 1,
  J = 10,
  N = 99,
  rho_null = 0,
  burnin = 500,
  rho_type = "pearson",
  del = 1e-10,
  trunc_lev = TRUE,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  errorType = "Gaussian",
  logNu = FALSE,
  sigvMethod = "factored",
  winsorize_eps = 0
)

Arguments

Value

An object of class "svp_test", a list containing:

s0

Test statistic from observed data.

sN

Simulated null distribution (vector of length N).

pval

Monte Carlo p-value.

test_type

Character string identifying the test.

null_param

Name of the parameter tested.

null_value

Value under the null hypothesis.

call

The matched call.

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2, leverage = TRUE, rho = -0.3)$y
test <- lmc_lev(y, p = 1, J = 10, N = 99)
print(test)

LMC Test for Student-t Tail Parameter

Description

Performs a Local Monte Carlo (LMC) test of the null hypothesis H_0: \nu = \nu_0 for the degrees of freedom parameter in an SV(p) model with Student-t errors. Testing \nu_0 = \infty (or a large value) corresponds to testing for normality.

Usage

lmc_t(
  y,
  p = 1,
  J = 10,
  N = 99,
  nu_null,
  burnin = 500,
  del = 1e-10,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  logNu = TRUE,
  direction = c("two-sided", "less", "greater"),
  sigvMethod = "factored",
  winsorize_eps = 0
)

Arguments

Value

An object of class "svp_test".

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2, errorType = "Student-t", nu = 5)$y
test <- lmc_t(y, p = 1, J = 10, N = 49, nu_null = 5)
print(test)

MMC Test for AR Order in SV(p) Models

Description

Performs a Maximized Monte Carlo (MMC) test of H_0: \phi_{p_0+1} = \cdots = \phi_p = 0 by maximizing the MC p-value over nuisance parameters (\phi_1, \ldots, \phi_{p_0}, \sigma_y, \sigma_v).

Usage

mmc_ar(
  y,
  p_null,
  p_alt,
  J = 10,
  N = 99,
  burnin = 500,
  eps = NULL,
  threshold = 1,
  method = "pso",
  maxit = NULL,
  del = 1e-10,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  errorType = "Gaussian",
  sigvMethod = "factored",
  logNu = TRUE,
  winsorize_eps = 0
)

Arguments

Value

A list with optimization output including value (maximized p-value) and par (nuisance parameters at the maximum).

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2)$y
mmc <- mmc_ar(y, p_null = 1, p_alt = 2, J = 10, N = 19,
              method = "pso", maxit = 10)
mmc$value

MMC Test for GED Shape Parameter

Description

Performs a Maximized Monte Carlo (MMC) test of H_0: \nu = \nu_0 for the GED shape parameter.

Usage

mmc_ged(
  y,
  p = 1,
  J = 10,
  N = 99,
  nu_null,
  burnin = 500,
  eps = NULL,
  threshold = 1,
  method = "pso",
  maxit = NULL,
  del = 1e-10,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  direction = c("two-sided", "less", "greater"),
  sigvMethod = "factored",
  winsorize_eps = 0
)

Arguments

Value

A list with optimization output including value (maximized p-value) and par (nuisance parameters at the maximum).

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2, errorType = "GED", nu = 1.5)$y
mmc <- mmc_ged(y, p = 1, J = 10, N = 19, nu_null = 2, method = "pso", maxit = 10)
mmc$value

MMC Test for Leverage in SV(p) Models

Description

Performs a Maximized Monte Carlo (MMC) test of the null hypothesis H_0: \rho = \rho_0 by maximizing the MC p-value over nuisance parameters (phi, sigma_y, sigma_v).

Usage

mmc_lev(
  y,
  p = 1,
  J = 10,
  N = 99,
  rho_null = 0,
  burnin = 500,
  eps = NULL,
  threshold = 1,
  method = "pso",
  maxit = NULL,
  rho_type = "pearson",
  del = 1e-10,
  trunc_lev = TRUE,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  errorType = "Gaussian",
  logNu = FALSE,
  sigvMethod = "factored",
  winsorize_eps = 0
)

Arguments

Value

A list with the optimization output including:

value

Maximized p-value.

par

Nuisance parameter values at the maximum.

Additional fields depend on the optimization method used.

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2, leverage = TRUE, rho = -0.3)$y
mmc <- mmc_lev(y, p = 1, J = 10, N = 19, method = "pso", maxit = 10)
mmc$value

MMC Test for Student-t Tail Parameter

Description

Performs a Maximized Monte Carlo (MMC) test of H_0: \nu = \nu_0 by maximizing the MC p-value over nuisance parameters (phi, sigma_y, sigma_v).

Usage

mmc_t(
  y,
  p = 1,
  J = 10,
  N = 99,
  nu_null,
  burnin = 500,
  eps = NULL,
  threshold = 1,
  method = "pso",
  maxit = NULL,
  del = 1e-10,
  wDecay = FALSE,
  Bartlett = FALSE,
  Amat = NULL,
  logNu = TRUE,
  direction = c("two-sided", "less", "greater"),
  sigvMethod = "factored",
  winsorize_eps = 0
)

Arguments

Value

A list with optimization output including value (maximized p-value) and par (nuisance parameters at the maximum).

Examples

y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2, errorType = "Student-t", nu = 5)$y
mmc <- mmc_t(y, p = 1, J = 10, N = 19, nu_null = 5, method = "pso", maxit = 10)
mmc$value

CDF of Standardized GED

Description

Computes the CDF of the standardized GED(\nu) distribution with unit variance.

Usage

pged_std(u, nu)

Arguments

Value

Numeric vector of probabilities.

Quantile function for standardized GED(nu) with Var = 1 Uses the relationship between GED CDF and incomplete gamma function.

Description

Quantile function for standardized GED(nu) with Var = 1 Uses the relationship between GED CDF and incomplete gamma function.

Usage

qged_std(p, nu)

Arguments

Value

Quantile value

Simulate from a Stochastic Volatility Model

Description

Master simulation function for SV(p) models. Supports Gaussian, Student-t, and GED error distributions, with optional leverage effects. This mirrors the interface of svp for estimation.

Usage

sim_svp(
  n,
  phi,
  sigy,
  sigv,
  errorType = "Gaussian",
  leverage = FALSE,
  rho = 0,
  nu = NULL,
  burnin = 500
)

Arguments

Details

The model is:

y_t = \sigma_y \exp(w_t / 2) z_t

w_t = \phi_1 w_{t-1} + \cdots + \phi_p w_{t-p} + \sigma_v v_t

where z_t follows a distribution specified by errorType (Gaussian, Student-t, or GED), and v_t is i.i.d. standard normal. When leverage = TRUE, the correlation between z_t and v_{t+1} is \rho.

For Student-t errors with leverage, the scale-mixture representation z_t = \zeta_t \lambda_t^{-1/2} is used, where leverage operates through the Gaussian component \zeta_t. For GED errors with leverage, a Gaussian copula construction z_t = F_{\mathrm{GED}}^{-1}(\Phi(\zeta_t)) is used. In both cases the returned z is the effective return innovation (not the latent \zeta_t), with marginal distribution matching the errorType.

Value

A named list of four length-n numeric vectors:

y

Observed returns y_t.

h

Log-volatility process w_t (equivalently h_t).

z

Return innovation such that y_t = \sigma_y \exp(h_t/2)\, z_t. Marginal distribution matches errorType: N(0,1) for Gaussian, t(\nu) for Student-t, unit-variance GED(\nu) for GED.

v

Volatility innovation such that h_t - \sum_{j=1}^p \phi_j h_{t-j} = \sigma_v\, v_t. Always N(0,1); under leverage, v_t = \rho\, \zeta_{t-1} + \sqrt{1-\rho^2}\, \epsilon_t.

See Also

svp for estimation.

Examples

# Gaussian SV(1), no leverage
sim <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2)
plot(sim$y, type = "l")

# Gaussian SV(1) with leverage
sim_lev <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2,
                   leverage = TRUE, rho = -0.5)
plot(sim_lev$y, type = "l")

# Student-t SV(1)
sim_t <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2,
                 errorType = "Student-t", nu = 5)
plot(sim_t$y, type = "l")

# GED SV(1)
sim_ged <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2,
                   errorType = "GED", nu = 1.5)
plot(sim_ged$y, type = "l")

Solve Discrete Lyapunov Equation

Description

Solves X = F X t(F) + Q for X using the vectorization approach.

Usage

solve_lyapunov_discrete(F_mat, Q)

Arguments

Value

Solution matrix X.

Estimate a Stochastic Volatility Model

Description

Master estimation function for SV(p) models using the Winsorized ARMA-SV (W-ARMA-SV) method. Supports Gaussian, Student-t, and GED error distributions, with optional leverage effects.

Usage

svp(
  y,
  p = 1,
  J = 10,
  leverage = FALSE,
  errorType = "Gaussian",
  rho_type = "pearson",
  del = 1e-10,
  trunc_lev = TRUE,
  wDecay = FALSE,
  logNu = FALSE,
  sigvMethod = "factored",
  winsorize_eps = 0
)

Arguments

Details

The model is:

y_t = \sigma_y \exp(w_t / 2) z_t

w_t = \phi_1 w_{t-1} + \cdots + \phi_p w_{t-p} + \sigma_v v_t

where z_t follows a distribution specified by errorType (Gaussian, Student-t, or GED), and v_t is i.i.d. standard normal. When leverage = TRUE, the correlation between z_t and v_t is estimated as \rho.

For Student-t errors with leverage, the correction factor C_t(\nu) from the scale-mixture representation is applied. For GED errors with leverage, the exact implicit equation is solved via 1D root-finding with Gauss-Hermite quadrature.

Value

Depending on errorType:

• "Gaussian": An object of class "svp" (see below).

• "Student-t": An object of class "svp_t".

• "GED": An object of class "svp_ged".

The "svp" class contains:

mu

Mean of \log(y_t^2 + \delta).

phi

Numeric vector of AR coefficients.

sigv

Standard deviation of volatility innovations.

sigy

Unconditional standard deviation.

rho

Leverage parameter (if estimated, otherwise NA).

y

The original data.

p, J

Model order and winsorizing parameter.

errorType

The error distribution used.

call

The matched call.

References

Ahsan, M. N. and Dufour, J.-M. (2021). Simple estimators and inference for higher-order stochastic volatility models. Journal of Econometrics, 224(1), 181-197. doi:10.1016/j.jeconom.2021.03.008

Ahsan, M. N., Dufour, J.-M., and Rodriguez-Rondon, G. (2026). Estimation and inference for stochastic volatility models with heavy-tailed distributions. Bank of Canada Staff Working Paper 2026-8. doi:10.34989/swp-2026-8

See Also

svpSE for standard errors.

Examples

# Gaussian SV(1) without leverage (default)
y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2)$y
fit <- svp(y)
summary(fit)

# With leverage
y2 <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2, leverage = TRUE, rho = -0.3)$y
fit2 <- svp(y2, leverage = TRUE)
coef(fit2)

# Student-t errors
y3 <- sim_svp(1000, phi = 0.9, sigy = 1, sigv = 0.2, errorType = "Student-t", nu = 5)$y
fit3 <- svp(y3, errorType = "Student-t")
summary(fit3)

Simulation-Based Standard Errors for SV(p) Models

Description

Computes standard errors and confidence intervals for estimated parameters by simulating from the fitted model and re-estimating. Supports all model types returned by svp: Gaussian (with or without leverage), Student-t, and GED.

Usage

svpSE(object, n_sim = 199, alpha = 0.05, burnin = 500, logNu = FALSE)

Arguments

Value

A list with:

CI

2 x k matrix of confidence intervals (lower, upper).

SEsim0

Standard errors relative to true parameter values.

SEsim

Standard errors relative to sample mean.

ISEconservative

Conservative interval-based standard errors.

ISEliberal

Liberal interval-based standard errors.

thetamat

Matrix of parameter estimates from simulations.

Examples

# Gaussian SV(1)
y <- sim_svp(1000, phi = 0.95, sigy = 1, sigv = 0.2)$y
fit <- svp(y)
se <- svpSE(fit, n_sim = 49)
se$CI

AR-order selection sweep for SV(p)

Description

Convenience wrapper around svp_IC: fits svp at each p = 1, ..., pmax and returns a matrix of information criteria along with the argmin per criterion.

Usage

svp_AR_order(
  y,
  pmax = 6L,
  J = 10L,
  leverage = FALSE,
  errorType = "Gaussian",
  rho_type = "pearson",
  del = 1e-10,
  trunc_lev = TRUE,
  wDecay = FALSE,
  logNu = FALSE,
  sigvMethod = "factored",
  winsorize_eps = 0L,
  filter_method = "mixture",
  proxy = c("bayes_optimal", "u"),
  K = 7L,
  M = 1000L,
  seed = 42L,
  criteria = c("BIC_Kalman", "AIC_Kalman", "BIC_HR", "AIC_HR")
)

Arguments

Value

A list with components:

IC

Numeric matrix, one row per criterion, one column per candidate p in 1:pmax.

argmin

Named integer vector, one entry per criterion, giving the selected p. NA_integer_ if all entries for that criterion are NA.

fits

List of length pmax containing the fitted svp() objects (or NULL if a fit failed).

See Also

svp_IC, svp, filter_svp

Examples

set.seed(1)
y <- sim_svp(2000, phi = 0.95, sigy = 1, sigv = 0.5)$y
res <- svp_AR_order(y, pmax = 4)
res$IC
res$argmin

Information criteria for SV(p) AR-order selection

Description

Computes information criteria for an svp fit to support AR-order selection. Eight criteria are computable; four are returned by default — BIC_Kalman, AIC_Kalman, BIC_HR, AIC_HR. These span two families (state-space QML and Hannan–Rissanen two-stage ARMA) and two penalty philosophies (Schwarz-consistent BIC / Shibata-efficient AIC), and were selected as the most informative criteria across the simulation grid of the SVHT methodology paper (Ahsan, Dufour and Rodriguez-Rondon 2026). The remaining four are available on request via the criteria argument.

Usage

svp_IC(
  fit,
  criteria = c("BIC_Kalman", "AIC_Kalman", "BIC_HR", "AIC_HR"),
  filter_method = c("mixture", "corrected", "particle"),
  proxy = c("bayes_optimal", "u"),
  K = 7L,
  M = 1000L,
  seed = 42L,
  del = 1e-10
)

Arguments

Details

Default criteria (returned by svp_IC(fit)):

• BIC_Kalman, AIC_Kalman: -2\,\hat\ell_K + k\log T and -2\,\hat\ell_K + 2k where \hat\ell_K is the (quasi-)log-likelihood from filter_svp; default filter_method = "mixture" uses the Gaussian mixture Kalman filter (Kim, Shephard and Chib 1998). BIC_Kalman is the primary recommended criterion: Schwarz-consistent under the Bayes-optimal leverage proxy (see proxy argument) and strong finite-sample performance across the simulation grid (Ahsan, Dufour and Rodriguez-Rondon 2026). AIC_Kalman is Shibata-efficient and often selects larger p sooner at p_{\mathrm{true}} \ge 2.

• BIC_HR, AIC_HR: Hannan–Rissanen (1982) two-stage ARMA(p,p) criteria. Stage 1: long-AR pre-whitening at order L = \lfloor 1.5\, T^{1/3}\rfloor produces residuals \hat\varepsilon_t. Stage 2: OLS regression of y_t^* on AR lags 1{:}p of y_t^* and MA lags 1{:}p of \hat\varepsilon_t gives \hat\sigma_u^2. Then T_{\mathrm{eff}} \log \hat\sigma_u^2 + \{2(2p{+}1), (2p{+}1)\log T_{\mathrm{eff}}\}. Filter-free anchor, robust to mis-specification of the GMKF mixture. BIC_HR is Schwarz-consistent for ARMA(p,p) (Hannan & Rissanen 1982; Pötscher 1989).

Opt-in criteria (request via criteria = ...):

• AICc_Kalman: AIC_Kalman with the Hurvich–Tsai (1989) small-sample correction 2k(k+1)/(T-k-1). Numerically equivalent to AIC_Kalman at T \ge 500; use when T < 500.

• BIC_Whittle: -2\,\hat\ell_W + k\log T where \hat\ell_W is the Whittle log-likelihood evaluated at the SV(p) signal-plus-noise spectral density f(\omega) = \sigma_v^2 / |1 - \sum_j \phi_j e^{-ij\omega}|^2 + \sigma_\varepsilon^2(\nu). Schwarz-consistent but collapses to \hat p = 1 in 98–100% of cells at p_{\mathrm{true}} \ge 2 under near-unit-root persistence (Ahsan, Dufour and Rodriguez-Rondon 2026). Useful as a conservative diagnostic: a Whittle selection of p > 1 is strong evidence against p = 1.

• AIC_YW, BIC_YW: Legacy / not recommended. Yule–Walker projection-error criteria on y_t^* = \log(y_t^2 + \delta) - \mu, computed as T \log \hat\sigma_{\mathrm{pred}}^2 + \{2k, k\log T\} with the AR(p) projection-error variance. Under SV(p), y_t^* is ARMA(p,p) (not AR(p)), so the AR projection error does not saturate at p_{\mathrm{true}} and the criteria are inconsistent: the AR(p) projection-error variance keeps decreasing past p_{\mathrm{true}}, producing non-monotone (sometimes anti-Schwarz) behaviour in T. Simulation evidence: 0–29% correct selection at p_{\mathrm{true}} = 2 across all DGP cells and T \le 10{,}000 (Ahsan, Dufour and Rodriguez-Rondon 2026). Retained for paper-reproducibility of the documented failure-case results; use BIC_HR / AIC_HR for theoretically consistent AR-order selection.

Lower is better; argmin over a grid of candidate p (see svp_AR_order) selects the AR order.

Value

Named numeric vector of the requested criteria. Lower is better.

Leverage invariance of non-Kalman criteria

Leverage does not affect AIC_YW, BIC_YW, or BIC_Whittle: under the W-ARMA-SV parameterization \mathrm{Cov}(v_t, \varepsilon_{t-1}) = 0 for all three error distributions (odd-times-even moment symmetry), so the autocovariance structure of y_t^* is invariant to the leverage parameter. The *_HR and *_Kalman criteria do incorporate leverage through the estimated \delta_p and the conditional state innovation variance.

References

Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control 19(6), 716–723. doi:10.1109/TAC.1974.1100705

Schwarz, G. (1978). Estimating the dimension of a model. Annals of Statistics 6(2), 461–464. doi:10.1214/aos/1176344136

Shibata, R. (1976). Selection of the order of an autoregressive model by Akaike's information criterion. Biometrika 63(1), 117–126. doi:10.1093/biomet/63.1.117

Hannan, E. J. (1980). The estimation of the order of an ARMA process. Annals of Statistics 8(5), 1071–1081. doi:10.1214/aos/1176345144

Hannan, E. J., and Rissanen, J. (1982). Recursive estimation of mixed autoregressive-moving average order. Biometrika 69(1), 81–94. doi:10.1093/biomet/69.1.81

Pötscher, B. M. (1989). Model selection under nonstationarity: Autoregressive models and stochastic linear regression models. Annals of Statistics 17(3), 1257–1274. doi:10.1214/aos/1176347267

Whittle, P. (1953). Estimation and information in stationary time series. Arkiv f\"or Matematik 2, 423–434. doi:10.1007/BF02590998

Dunsmuir, W. (1979). A central limit theorem for parameter estimation in stationary vector time series and its application to models for a signal observed with noise. Annals of Statistics 7(3), 490–506. doi:10.1214/aos/1176344671

Hurvich, C. M., and Tsai, C.-L. (1989). Regression and time series model selection in small samples. Biometrika 76(2), 297–307. doi:10.1093/biomet/76.2.297

Kim, S., Shephard, N., and Chib, S. (1998). Stochastic volatility: likelihood inference and comparison with ARCH models. Review of Economic Studies 65(3), 361–393. doi:10.1111/1467-937X.00050

White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica 50(1), 1–25. doi:10.2307/1912526

Ahsan, M. N., Dufour, J.-M., and Rodriguez-Rondon, G. (2026). Estimation and inference for stochastic volatility models with heavy-tailed distributions. Bank of Canada Staff Working Paper 2026-8. doi:10.34989/swp-2026-8

See Also

svp_AR_order, svp, filter_svp

Examples

set.seed(1)
y <- sim_svp(2000, phi = 0.95, sigy = 1, sigv = 0.5)$y
fit1 <- svp(y, p = 1)
fit2 <- svp(y, p = 2)
svp_IC(fit1)
svp_IC(fit2)