sstvars

install.packages("sstvars")

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

# 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)