Learn R Programming
⚠️There's a newer version (1.1.6) of this package.Take me there.

sstvars

The goal of sstvars is to provide a comprehensive toolkit for maximum likelihood (ML) estimation and analysis of reduced form and structural smooth transition vector autoregressive (STVAR) models (including threshold VAR models). Various transition weight functions, conditional distributions, and identification methods are accommodated. Also constrained ML estimation is supported with constraints on the autoregressive parameters, regimewise means, weight parameters, and the impact matrix. See the vignette for a more detailed description of the package.

Installation

You can install the released version of gmvarkit from CRAN with:

install.packages("sstvars")

You can install the development version of sstvars from GitHub with:

# install.packages("devtools")
devtools::install_github("saviviro/sstvars")

Example

This is a basic example on how to use sstvars in time series analysis. The estimation process is computationally demanding and takes advantage of parallel computing. After estimating the model, it is shown by simple examples how to conduct some further analysis.

# These examples use the data 'gdpdef' which comes with the package, and contains the quarterly percentage growth rate
# of real U.S. GDP and quarterly percentage growth rate of U.S. GDP implicit price deflator, covering the period 
# from 1959Q1 to 2019Q4.
data(gdpdef, package="sstvars")

# Some of the below examples are computationally demanding. Running them all will take approximately 10 minutes.

### Reduced form STVAR models ###

# Estimate a reduced form two-regime Student's t STVAR p=2 model with threshold transition weight function using the first
# lag of the first variable (GDP) as the switching variable. The below estimation is based on two estimation
# rounds with seeds set for reproducibility.
# (IMPORTANT: typically empirical applications require more estimation rounds, e.g., tens, hundreds or even thousand, depending
# on the size of the model, and with the two-phase procedure often much more).
fit <- fitSTVAR(gdpdef, p=2, M=2, weight_function="threshold", weightfun_pars=c(2, 1), cond_dist="Student",
                estim_method="three-phase", nrounds=2, ncores=2, seeds=1:2)
                
# Information on the estimated model:
plot(fit) # Plot the estimated transition weight function with data
summary(fit) # Summary printoout of the estimated model
get_foc(fit) # The first order condition (gradient of the log-likelihood function)
get_soc(fit) # The second order condition (eigenvalues of approximated Hessian)
profile_logliks(fit) # Plot profile log-likelihood functions about the estimate

# See also the functions alt_stvar and filter_estimates.

# Check the stationarity condition for the estimated model, i.e., that the 
# upper bound of the joint spectral radius is less than one:
bound_JSR(fit, epsilon=0.1, ncores=2) # Adjust epsilon for a tighter bound

# Estimate the above model but with the autoregressive matrices restricted to be equal in both regimes
# (so that only the intercepts and the conditional covariance matrix vary in time):
C_mat <- rbind(diag(2*2^2), diag(2*2^2))
fitc <- fitSTVAR(gdpdef, p=2, M=2, weight_function="threshold", weightfun_pars=c(2, 1), cond_dist="Student",
                 AR_constraints=C_mat, nrounds=2, ncores=2, seeds=1:2)

# Estimate the above model but with the autoregressive matrices and unconditional means restricted to be equal
# in both regimes (so that only the conditional covariance matrix varies in time), two-phase estimation 
# is used because mean-constraints are not supported in the three-phase estimation:
fitcm <- fitSTVAR(gdpdef, p=2, M=2, weight_function="threshold", weightfun_pars=c(2, 1), cond_dist="Student",
                  AR_constraints=C_mat, mean_constraints=list(1:2), nrounds=2, ncores=2, seeds=1:2)

# Estimate the above logistic STVAR model without constraints on the autoregressive parameters but with the 
# the location parameter constrained to 1 and scale parameter unconstrained. Two-phase estimation is used because
# this type of weight constraints  are not supported in the three-phase estimation:
fitw <- fitSTVAR(gdpdef, p=2, M=2, weight_function="logistic", weightfun_pars=c(2, 1), cond_dist="Student",
                 weight_constraints=list(R=matrix(c(0, 1), nrow=2), r=c(1, 0)), nrounds=2, ncores=2, seeds=1:2)

# The constraints can be tested with the functions LR_test, Wald_test, and Rao_test.

# Residual based model diagnostics:
diagnostic_plot(fit, type="series", resid_type="standardized") # Standardized residual time series
diagnostic_plot(fit, type="ac", resid_type="raw") # Autocorrelation function of unstandardized residuals
diagnostic_plot(fit, type="ch", resid_type="standardized") # Autocorrelation function of squared standardized residuals
diagnostic_plot(fit, type="dist", resid_type="standardized") # Histograms and Q-Q plots of standardized residuals

Portmanteau_test(fit, nlags=20, which_test="autocorr") # Portmanteau test for remaining autocorrelation
Portmanteau_test(fit, nlags=20, which_test="het.sked") # Portmanteau test applied for testing cond. het.kedasticity

# Simulate a sample path from the estimated model, initial drawn from the first regime:
sim <- simulate(fit, nsim=100, init_regime=1)

# Forecast future values of the process:
pred <- predict(fit, nsteps=10, ci=c(0.95, 0.80))
plot(pred)


### Structural STVAR models ###

# stvars implements two identification methods: recursive identification and
# identification by heteroskedasticity. The structural models are estimated 
# based on preliminary estimates from a reduced form model. If the structural model
# is not overidentifying, the model is merely reparametrized and no estimation is
# required (recursively identified models are just reduced form models marked as structural). 

# Identify the above threshold VAR model by recursive identification:
fitrec <- fitSSTVAR(fit, identification="recursive")
fitrec

# Identify the above threshold VAR model by heteroskedasticity:
fithet <- fitSSTVAR(fit, identification="heteroskedasticity")
fithet

# Identification by non-Gaussianity available for models with independent Student's t distribution
# or independent skewed t distribution as the conditional distribution. The reduced form model is
# then readily identified by non-Gaussianity. Estimate a reduced form model identified by
# non-Gaussianity with independent Student's t shocks: 
fitindt <- fitSTVAR(gdpdef, p=1, M=2, weight_function="logistic", weightfun_pars=c(1, 1),
                    cond_dist="ind_Student", nrounds=2, ncores=2, seeds=1:2)
fitindt

# Impose overidentying constraint with the argument B_constraints by estimating
# with fitSSTVARs:
fitindtb <- fitSSTVAR(fitindt, identification="non-Gaussianity",
                      B_constraints=matrix(c(NA, NA, 0, NA), nrow=2))

# Reorder the columns of the impact matrix of fithet to the reverse ordering:
fithet <- reorder_B_columns(fithet, perm=c(2, 1))
fithet

# Change all signs of the first column of the impact matrix of fithet:
fithet <- swap_B_signs(fithet, which_to_swap=1)
fithet

# Structural models based on different orderings of signs of the columns of any
# single B_1,...B_M can be estimated by specifying the arguments B_pm_reg, B_perm,
# and B_signs in fitSSTVAR. 


# Estimate the generalized impulse response function (GIRF) for the recursively
# identified model to one-standard-error positive shocks with the starting values
# generated form the first regime, N=30 steps ahead and 95% confidence intervals 
# that reflect uncertainty about the initial value within the regime:
girf1 <- GIRF(fitrec, which_shocks=1:2, shock_size=1, N=30, init_regime=1, ci=0.95)
plot(girf1)

# Estimate the above GIRF but instead of drawing initial values form the first regime,
# use tha last p observations of the data as the initial values:
girf2 <- GIRF(fitrec, which_shocks=1:2, shock_size=1, N=30, init_values=fitrec$data)
plot(girf2)

# Estimate the generalized impulse response function (GIRF) for the recursively
# identified model to two-standard-error negative shocks with the starting values
# generated form the second regime, N=30 steps ahead and 95% confidence intervals 
# that reflect uncertainty about the initial value within the regime. Also, scale
# the responses to the first shock to correspond to a 0.3 increase of the first variable.
# Moreover, accumulate the responses of the second variable.
girf3 <- GIRF(fitrec, which_shocks=1:2, shock_size=-2, N=30, init_regime=2, 
              scale=c(1, 1, 0.3), which_cumulative=2, ci=0.95)
plot(girf3)

# Estimate the generalized forecast error variance decomposition (GFEVD) for the 
# recursively identified model with the initial values being all possible length p
# histories in the data, N=30 steps ahead to one-standard-error positive shocks. 
gfevd1 <- GFEVD(fitrec, shock_size=1, N=30, initval_type="data")
plot(gfevd1)

# Estimate the GFEVD for the recursively identified model with the initial values
# being all possible length p histories in the data AND the signs and sizes of the
# corresponding shocks being the identified structural shocks recovered from the
# fitted model.
gfevd2 <- GFEVD(fitrec, N=30, use_data_shocks=TRUE)
plot(gfevd2)

# Plot time series of the impact response GFEVDs to the first shock to examine 
# the contribution of the first shocks to the forecast error variances at impact
# at each point of time:
plot(gfevd2, data_shock_pars=c(1, 0))

# Estimate the linear impulse response function for the recursively identified
# STVAR model based on the first regime, the responses of the second variable
# accumulated:
irf1 <- linear_IRF(fitrec, N=30, regime=1, which_cumulative=2)
plot(irf1)
# The above model is nonlinear, so confidence bounds (that reflect the uncertainty
# about the parameter estimate) are not available.

# Bootstrapped confidence bounds can be calculated for models with time-invariant
# autoregressive coeffients, e.g., the restricted model fitcm estimated above. 
# Identify the shocks if fitcm by heteroskedasticity:
fitcmhet <- fitSSTVAR(fitcm, identification="heteroskedasticity")
fitcmhet

# Estimate the linear impulse reponse function for fitcmhet with bootstrapped
# 95% confidence bounds that reflect uncertainty about the parameter estimates:
irf2 <- linear_IRF(fitcmhet, N=30, ci=0.95, bootstrap_reps=250)
plot(irf2)

References

  • Anderson H., Vahid F. (1998) Testing multiple equation systems for common nonlinear components. Journal of Econometrics, 84:1, 1-36.
  • Hubrich K., Teräsvirta. T. (2013). Thresholds and Smooth Transitions in Vector Autoregressive Models. CREATES Research Paper 2013-18, Aarhus University.
  • Kheifets I., Saikkonen P. (2020). Stationarity and ergodicity of vector STAR models. Econometrics Review, 39:407-414, 1311-1324.
  • Koop G., Pesaran M.H., Potter S.M. (1996). Impulse response analysis in nonlinear multivariate models. Journal of Econometrics, 74:1, 119-147.
  • Lanne M., Virolainen S. 2024. A Gaussian smooth transition vector autoregressive model: An application to the macroeconomic effects of severe weather shocks. Unpublished working paper, available as arXiv:2403.14216.
  • Virolainen S. 2024. Identification by non-Gaussianity in structural threshold and smooth transition vector autoregressive. Unpublished working paper, available as arXiv:2404.19707.

Copy Link

Version

Install

install.packages('sstvars')

Monthly Downloads

440

Version

1.1.0

License

GPL-3

Issues

1

Pull Requests

0

Stars

5

Forks

2

Maintainer

Savi Virolainen

Last Published

November 29th, 2024

Functions in sstvars (1.1.0)

Rao_test

Perform Rao's score test for a STVAR model
Student_densities_Cpp

Calculate log multivariate Student's t densities
Portmanteau_test

Perform adjusted Portmanteau test for a STVAR model
GAfit

Genetic algorithm for preliminary estimation of reduced form STVAR models
STVAR

Create a class 'stvar' object defining a reduced form or structural smooth transition VAR model
GFEVD

Estimate generalized forecast error variance decomposition for structural STVAR models.
GIRF

Estimate generalized impulse response function for structural STVAR models.
LR_test

Perform likelihood ratio test for a STVAR model
Gaussian_densities_Cpp

Calculate log multivariate Gaussian densities
Gaussian_densities_const_Cpp

Calculate log multivariate Gaussian densities
Wvec

Vectorization operator that removes zeros
alt_stvar

Construct a STVAR model based on results from an arbitrary estimation round of fitSTVAR
acidata

A monthly U.S. data covering the period from 1961I to 2022III (735 observations) and consisting four variables. First, The Actuaries Climate Index (ACI), which is a measure of the frequency of severe weather and the extend changes in sea levels. Second, the monthly GDP growth rate constructed by the Federal Reserve Bank of Chicago from a collapsed dynamic factor analysis of a panel of 500 monthly measures of real economic activity and quarterly real GDP growth. Third, the monthly growth rate of the consumer price index (CPI). Third, an interest rate variable, which is the effective federal funds rate that is replaced by the the Wu and Xia (2016) shadow rate during zero-lower-bound periods. The Wu and Xia (2016) shadow rate is not bounded by the zero lower bound and also quantifies unconventional monetary policy measures, while it closely follows the federal funds rate when the zero lower bound does not bind.
all_pos_ints

Check whether all arguments are positive integers
calc_gradient

Calculate gradient or Hessian matrix
bounding_const_M

Compute the bounding constant for acceptance-rejection sampling
Wald_test

Perform Wald test for a STVAR model
VAR_pcovmat

Calculate the dp-dimensional covariance matrix of p consecutive observations of a VAR process
bound_jsr_G

Calculate upper bound for the joint spectral radius of a set of matrices
bound_JSR

Calculate upper bound for the joint spectral radius of the "companion form AR matrices" of the regimes
check_weightfun_pars

Check the argument weightfun_pars
check_Bt_Cpp

Check Matrix B Invertibility with C++ (Internal Function)
check_data

Check the data is in the correct form
check_constraints

Check the constraint matrix has the correct form
check_stvar

Checks whether the given object has class attribute 'stvar'
check_pMd

Check that p, M, and d are correctly set
change_regime

Change the parameters of a specific regime of the given parameter vector
change_parametrization

Change parametrization of a parameter vector
check_params

Check whether the parameter vector is in the parameter space and throw error if not
check_exoweights

Checks whether the given exogenous transition weights for simulation are correctly specified.
diag_Omegas

Simultaneously diagonalize two covariance matrices
diagnostic_plot

Residual diagnostic plot for a STVAR model
create_J_matrix

Create a special matrix J
filter_estimates

Filter inappropriate the estimates produced by fitSTVAR
fitSSTVAR

Maximum likelihood estimation of a structural STVAR model based on preliminary estimates from a reduced form model.
fitSTVAR

Two-phase or three-phase (penalized) maximum likelihood estimation of a reduced form smooth transition VAR model
create_Fi_matrix

Create Matrix F_i
estim_LS

Internal estimation function for estimating autoregressive and threshold parameters of TVAR models by the method of least squares.
fitbsSSTVAR

Internal estimation function for estimating STVAR model when bootstrapping confidence bounds for IRFs in linear_IRF
estim_NLS

Internal estimation function for estimating autoregressive and weight parameters of STVAR models by the method of nonlinear least squares.
format_valuef

Function factory for value formatting
get_boldA_eigens

Calculate absolute values of the eigenvalues of the "bold A" matrices containing the AR coefficients
get_IC

Calculate AIC, HQIC, and BIC
get_alpha_mt

Get the transition weights alpha_mt
generate_skewed_t

Generate random samples from the skewed t-distribution
get_boldA_eigens_par

Calculate absolute values of the eigenvalues of the "bold A" matrices containing the AR coefficients
gdpdef

U.S. real GDP percent change and GDP implicit price deflator percent change.
form_boldA

Form the \((dp\times dp)\) "bold A" matrices related to the VAR processes
get_Bt_Cpp

Calculate the impact matrix \(B_t\) for all \(t\) for models with a non-Gaussian conditional distribution with mutually independent shocks.
get_Sigmas

Calculate the dp-dimensional covariance matrices \(\Sigma_{m,p}\) in the transition weights with weight_function="relative_dens"
get_hetsked_sstvar

Switch from two-regime reduced form STVAR model to a structural model identified by heteroskedasticity
get_residuals

Calculate residuals of a smooth transition VAR
get_mu_yt_Cpp

Calculate the conditional means of the process
get_omega_eigens

Calculate the eigenvalues of the "Omega" error term covariance matrices
get_omega_eigens_par

Calculate the eigenvalues of the "Omega" error term covariance matrices
get_symmetric_sqrt

Calculate symmetric square root matrix of a positive definite covariance matrix
get_regime_autocovs

Calculate regimewise autocovariance matrices
get_new_start

Get the new starting time of series that is forwarded some number of steps
get_minval

Returns the default smallest allowed log-likelihood for given data.
get_regime_means

Calculate regime means \(\mu_{m}\)
linear_IRF

Estimate linear impulse response function based on a single regime of a structural STVAR model.
mat_power

Compute the j:th power of a square matrix A
in_paramspace

Determine whether the parameter vector is in the parameter space
iterate_more

Maximum likelihood estimation of a reduced form or structural STVAR model based on preliminary estimates
order_B

Reorder columns of a square matrix so that the first nonzero elements are in decreasing order
pick_Am

Pick coefficient matrices
ind_Student_densities_Cpp

Calculate log independent multivariate Student's t densities
ind_skewed_t_densities_Cpp

Calculate log independent multivariate skewed t densities
n_params

Calculate the number of (freely estimaed) parameters in the model
loglikelihood

Log-likelihood function
pick_W

Pick the structural parameter matrix W
plot_struct_shocks

Plot structural shock time series of a STVAR model
pick_phi0

Pick \(\phi_{m,0}\) or \(\mu_{m}\), m=1,..,M vectors
pick_Omegas

Pick covariance matrices
pick_Ami

Pick coefficient matrix
pick_lambdas

Pick the structural parameter eigenvalues 'lambdas'
pick_distpars

Pick distribution parameters
pick_allA

Pick all coefficient matrices
pick_regime

Pick regime parameters
pick_weightpars

Pick transition weight parameters
print.hypotest

Print method for the class hypotest
profile_logliks

Plot profile log-likelihood functions about the estimates
print.stvarsum

Summary print method from objects of class 'stvarsum'
random_coefmats2

Create random stationary VAR model \((dxd)\) coefficient matrices \(A\).
random_covmat

Create random VAR model error term covariance matrix
random_distpars

Create random distribution parameter values
random_coefmats

Create random VAR model \((dxd)\) coefficient matrices \(A\).
plot.stvarpred

Predict method for class 'stvar' objects
random_impactmat

Create random VAR model impact matrix
random_ind

Create random mean parametrized parameter vector
random_weightpars

Create random transition weight parameter values
reform_constrained_pars

Reform constrained parameter vector into the "standard" form
regime_distance

Calculate "distance" between two (scaled) regimes \(\upsilon_{m}\)\( = (\phi_{m,0},\)\(\phi_{m}\)\(,\sigma_{m})\)
skewed_t_dens

The density function of the univariate skewed t distribution
smart_covmat

Create random VAR model \((dxd)\) error term covariance matrix \(\Omega\) fairly close to the given positive definite covariance matrix using (scaled) Wishart distribution
reorder_B_columns

Reorder columns of impact matrix B (and lambda parameters if any) of a structural STVAR model that is identified by heteroskedasticity or non-Gaussianity.
reform_data

Reform data
simulate.stvar

Simulate method for class 'stvar' objects
simulate_from_regime

Simulate observations from a regime of a STVAR model
redecompose_Omegas

In the decomposition of the covariance matrices (Muirhead, 1982, Theorem A9.9), change the ordering of the covariance matrices.
sort_impactmats

Sort and sign change the columns of the impact matrices of the regimes so that the first element in each column of \(B_1\) is positive and in a decreasing order.
sstvars-package

sstvars: toolkit for reduced form and structural smooth transition vector autoregressive models
smart_weightpars

Create random transition weight parameter values
sort_regimes

Sort regimes in parameter vector according to transition weights into a decreasing order
smart_ind

Create random parameter vector that is fairly close to a given parameter vector
stab_conds_satisfied

Check the stability condition for each of the regimes
smart_distpars

Create random distribution parameter values close to given values
smart_impactmat

Create a random VAR model \((dxd)\) error impact matrix \(B\) fairly close to the given invertible impact matrix.
swap_parametrization

Swap the parametrization of a STVAR model
swap_B_signs

Swap all signs in pointed columns of the impact matrix of a structural STVAR model that is identified by heteroskedasticity or non-Gaussianity
stvar_to_sstvars110

Update STVAR model estimated with a version of the package <1.1.0 to be compatible with the versions >=1.1.0.
standard_errors

Calculate standard errors for estimates of a smooth transition VAR model
vech

Parsimonious vectorization operator for symmetric matrices
warn_eigens

Warn about near-unit-roots in some regimes
stand_t_dens

The density function of the univariate t distribution with zero mean and unit variance
unWvec

Reverse vectorization operator that restores zeros
uncond_moments

Calculate the unconditional means, variances, the first p autocovariances, and the first p autocorrelations of the regimes of the model.
unvec

Reverse vectorization operator
usamone

A quarterly U.S. data covering the period from 1954Q3 to 2021Q4 (270 observations) and consisting three variables: cyclical component of the log of real GDP, the log-difference of GDP implicit price deflator, and an interest rate variable. The interest rate variable is the effective federal funds rate from 1954Q3 to 2008Q2 and after that the Wu and Xia (2016) shadow rate, which is not constrained by the zero lower bound and also quantifies unconventional monetary policy measures. The log-differences of the GDP deflator and producer price index are multiplied by hundred.
vec

Vectorization operator
unvech

Reverse operator of the parsimonious vectorization operator vech
usacpu

A monthly U.S. data covering the period from 1987:4 to 2024:2 (443 observations) and consisting six variables. First, the climate policy uncertainty index (CPUI) (Gavridiilis, 2021), which is a news based measure of climate policy uncertainty. Second, the economic policy uncertainty index (EPUI), which is a news based measure of economic policy uncertainty. Third, the log-difference of real indsitrial production index (IPI). Fourth, the log-difference of the consumer price index (CPI). Fifth, the log-difference of the producer price index (PPI). Sixth, an interest rate variable, which is the effective federal funds rate that is replaced by the the Wu and Xia (2016) shadow rate during zero-lower-bound periods. The Wu and Xia (2016) shadow rate is not bounded by the zero lower bound and also quantifies unconventional monetary policy measures, while it closely follows the federal funds rate when the zero lower bound does not bind. This is the dataset used in Virolainen (2024)