Package 'sstvars'

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

Description

sstvars is a package for reduced form and structural smooth transition vector autoregressive models. The package implements various transition weight functions, conditional distributions, identification methods, and parameter restrictions. The model parameters are estimated with the method of maximum likelihood by running multiple rounds of a two-phase estimation procedure in which a genetic algorithm is used to find starting values for a gradient based method. For evaluating the adequacy of the estimated models, sstvars utilizes residuals based diagnostics and provides functions for graphical diagnostics and for calculating formal diagnostic tests. sstvars also accommodates the estimation of linear impulse response functions, nonlinear generalized impulse response functions, and generalized forecast error variance decompositions. Further functionality includes hypothesis testing, plotting the profile log-likelihood functions about the estimate, simulation from STVAR processes, and forecasting, for example.

The vignette is a good place to start, and see also the readme file.

Author(s)

you <[email protected]>

See Also

Useful links:

• https://github.com/saviviro/sstvars

• Report bugs at https://github.com/saviviro/sstvars/issues

---

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.

Description

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.

Usage

acidata

Format

A numeric matrix of class 'ts' with 735 rows and 4 columns with one time series in each column:

First column (GDP):

The cyclical component of the log of real GDP, https://fred.stlouisfed.org/series/GDPC1.

Second column (GDPDEF):

The log-difference of GDP implicit price deflator, https://fred.stlouisfed.org/series/GDPDEF.

Third column (RATE):

The Federal funds rate from 1954Q3 to 2008Q2 and after that the Wu and Xia (2016) shadow rate, https://fred.stlouisfed.org/series/FEDFUNDS, https://www.atlantafed.org/cqer/research/wu-xia-shadow-federal-funds-rate.

Source

The Federal Reserve Bank of St. Louis database and the Federal Reserve Bank of Atlanta's website

References

• American Academy of Actuaries, Canadian Institute of Actuaries, Casualty Actuarial Society, and Society of Actuaries, 2023. Actuaries Climate Index. https://actuariesclimateindex.org.

• Federal Reserve Bank of Chicago, 2023. Monthly GDP Growth Rate Data. https://www.chicagofed.org/publications/bbki/index.

• Wu J. and Xia F. 2016. Measuring the macroeconomic impact of monetary policy at the zero lower bound. Journal of Money, Credit and Banking, 48(2-3): 253-291.

---

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

Description

alt_stvar constructs a STVAR model based on results from an arbitrary estimation round of fitSTVAR

Usage

alt_stvar(stvar, which_largest = 1, which_round, calc_std_errors = FALSE)

Arguments

stvar	object of class "stvar"
which_largest	based on estimation round with which largest log-likelihood should the model be constructed? An integer value in 1,...,nrounds. For example, which_largest=2 would take the second largest log-likelihood and construct the model based on the corresponding estimates.
which_round	based on which estimation round should the model be constructed? An integer value in 1,...,nrounds. If specified, then which_largest is ignored.
calc_std_errors	should approximate standard errors be calculated?

Details

It's sometimes useful to examine other estimates than the one with the highest log-likelihood. This function is wrapper around STVAR that picks the correct estimates from an object returned by fitSTVAR.

Value

Returns an S3 object of class 'stvar' defining a smooth transition VAR model. The returned list contains the following components (some of which may be NULL depending on the use case):

data	The input time series data.
model	A list describing the model structure.
params	The parameters of the model.
std_errors	Approximate standard errors of the parameters, if calculated.
transition_weights	The transition weights of the model.
regime_cmeans	Conditional means of the regimes, if data is provided.
total_cmeans	Total conditional means of the model, if data is provided.
total_ccovs	Total conditional covariances of the model, if data is provided.
uncond_moments	A list of unconditional moments including regime autocovariances, variances, and means.
residuals_raw	Raw residuals, if data is provided.
residuals_std	Standardized residuals, if data is provided.
structural_shocks	Recovered structural shocks, if applicable.
loglik	Log-likelihood of the model, if data is provided.
IC	The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.
all_estimates	The parameter estimates from all estimation rounds, if applicable.
all_logliks	The log-likelihood of the estimates from all estimation rounds, if applicable.
which_converged	Indicators of which estimation rounds converged, if applicable.
which_round	Indicators of which round of optimization each estimate belongs to, if applicable.
LS_estimates	The least squares estimates of the parameters in the form $(\phi_{1,0},...,\phi_{M,0},\varphi_1,...,\varphi_M,\alpha$ (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

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.

• 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.

• Kheifets I.L., Saikkonen P.J. 2020. Stationarity and ergodicity of Vector STAR models. Econometric Reviews, 39:4, 407-414.

• Lütkepohl H., Netšunajev A. 2017. Structural vector autoregressions with smooth transition in variances. Journal of Economic Dynamics & Control, 84, 43-57.

• Tsay R. 1998. Testing and Modeling Multivariate Threshold Models. Journal of the American Statistical Association, 93:443, 1188-1202.

• Virolainen S. 2024. Identification by non-Gaussianity in structural threshold and smooth transition vector autoregressive models. Unpublished working paper, available as arXiv:2404.19707.

See Also

STVAR

Examples

## These are long-running examples that take approximately 10 seconds to run.

# Estimate a Gaussian STVAR p=1, M=2 model with exponential weight function and
# the first lag of the second variable as the switching variables. Run only two
# estimation rounds:
fit12 <- fitSTVAR(gdpdef, p=1, M=2, weight_function="exponential", weightfun_pars=c(2, 1),
 nrounds=2, seeds=c(1, 7))
fit12$loglik # Log-likelihood of the estimated model

# Print the log-likelihood obtained from each estimation round:
fit12$all_logliks

# Construct the model based on the second largest log-likelihood found in the
# estimation procedure:
fit12_alt <- alt_stvar(fit12, which_largest=2, calc_std_errors=FALSE)
fit12_alt$loglik # Log-likelihood of the alternative solution

# Construct a model based on a specific estimation round, the second round:
fit12_alt2 <- alt_stvar(fit12, which_round=2, calc_std_errors=FALSE)
fit12_alt2$loglik # Log-likelihood of the alternative solution

---

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

Description

bound_JSR calculates an bounds for the joint spectral radius of the "companion form AR matrices" matrices of the regimes to assess the validity of the stationarity condition.

Usage

bound_JSR(
  stvar,
  epsilon = 0.02,
  adaptive_eps = FALSE,
  ncores = 2,
  print_progress = TRUE
)

Arguments

stvar	object of class "stvar"
epsilon	a strictly positive real number that approximately defines the goal of length of the interval between the lower and upper bounds. A smaller epsilon value results in a narrower interval, thus providing better accuracy for the bounds, but at the cost of increased computational effort. Note that the bounds are always wider than epsilon and it is not obvious what epsilon should be chosen obtain bounds of specific tightness.
adaptive_eps	logical: if TRUE, starts with a large epsilon and then decreases it gradually whenever the progress of the algorithm requires, until the value given in the argument epsilon is reached. Usually speeds up the algorithm substantially but is an unconventional approach, and there is no guarantee that the algorithm converges appropriately towards bounds with the tightness given by the argument epsilon.
ncores	the number of cores to be used in parallel computing.
print_progress	logical: should the progress of the algorithm be printed?

Details

A sufficient condition for ergodic stationarity of the STVAR processes implemented in sstvars is that the joint spectral radius of the "companion form AR matrices" of the regimes is smaller than one (Kheifets and Saikkonen, 2020). This function calculates an upper (and lower) bound for the JSR and is implemented to assess the validity of this condition in practice. If the bound is smaller than one, the model is deemed ergodic stationary.

Implements the branch-and-bound method by Gripenberg (1996) in the conventional form (adaptive_eps=FALSE) and in a form incorporating "adaptive tightness" (adaptive_eps=FALSE). The latter approach is unconventional and does not guarantee appropriate convergence of the bounds close to the desired tightness given in the argument epsilon, but it usually substantially speeds up the algorithm. When print_progress==TRUE, the tightest bounds found so-far are printed in each iteration of the algorithm, so you can also just terminate the algorithm when the bounds are tight enough for your purposes. Consider also adjusting the argument epsilon, in particular when adaptive_eps=FALSE, as larger epsilon does not just make the bounds less tight but also speeds up the algorithm significantly. See Chang and Blondel (2013) for a discussion on variuous methods for bounding the JSR.

Value

Returns lower and upper bounds for the joint spectral radius of the "companion form AR matrices" of the regimes.

References

• C-T Chang and V.D. Blondel. 2013 . An experimental study of approximation algorithms for the joint spectral radius. Numerical algorithms, 64, 181-202.

• Gripenberg, G. 1996. Computing the joint spectral radius. Linear Algebra and its Applications, 234, 43–60.

• I.L. Kheifets, P.J. Saikkonen. 2020. Stationarity and ergodicity of Vector STAR models. Econometric Reviews, 39:4, 407-414.

• Virolainen S. 2024. Identification by non-Gaussianity in structural threshold and smooth transition vector autoregressive models. Unpublished working paper, available in ArXiv.

See Also

bound_jsr_G

Examples

## Below examples take approximately 5 seconds to run.

# Gaussian STVAR p=1, M=2 model with weighted relative stationary densities
# of the regimes as the transition weight function:
theta_122relg <- c(0.734054, 0.225598, 0.705744, 0.187897, 0.259626, -0.000863,
  -0.3124, 0.505251, 0.298483, 0.030096, -0.176925, 0.838898, 0.310863, 0.007512,
  0.018244, 0.949533, -0.016941, 0.121403, 0.573269)
mod122 <- STVAR(data=gdpdef, p=1, M=2, params=theta_122relg)

# Absolute values of the eigenvalues of the "companion form AR matrices":
summary(mod122)$abs_boldA_eigens
# It is a necessary (but not sufficient!) condition for ergodic stationary that
# the spectral radius of the "companion form AR matrices" are smaller than one
# for all of the regimes. A sufficient (but not necessary) condition for
# ergodic stationary is that the joint spectral radius of the companion form
# AR matrices" of the regimes is smaller than one. Therefore, we calculate
# bounds for the joint spectral radius.

## Bounds by Gripenberg's (1996) branch-and-bound method:
# Since the largest modulus of the companion form AR matrices is not very close
# to one, we likely won't need very thight bounds to verify the JSR is smaller
# than one. Thus, using a small epsilon would make the algorithm unnecessarily slow,
# so we use the (still quite small) epsilon=0.01:
bound_JSR(mod122, epsilon=0.01, adaptive_eps=FALSE)
# The upper bound is smaller than one, so the model is ergodic stationary.

# If we want tighter bounds, we can set smaller epsilon, e.g., epsilon=0.001:
bound_JSR(mod122, epsilon=0.001, adaptive_eps=FALSE)

# Using adaptive_eps=TRUE usually speeds up the algorithm when the model
# is large, but with the small model here, the speed-difference is small:
bound_JSR(mod122, epsilon=0.001, adaptive_eps=TRUE)

---

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

Description

bound_jsr_G calculates lower and upper bounds for the joint spectral radious of a set of square matrices, typically the "bold A" matrices, using the algorithm by Gripenberg (1996).

Usage

bound_jsr_G(
  S,
  epsilon = 0.01,
  adaptive_eps = FALSE,
  ncores = 2,
  print_progress = TRUE
)

Arguments

S	the set of matrices the bounds should be calculated for in an array, in STVAR applications, all $((dp)x(dp))$ "bold A" (companion form) matrices in a 3D array, so that [, , m] gives the matrix the regime m.
epsilon	a strictly positive real number that approximately defines the goal of length of the interval between the lower and upper bounds. A smaller epsilon value results in a narrower interval, thus providing better accuracy for the bounds, but at the cost of increased computational effort. Note that the bounds are always wider than epsilon and it is not obvious what epsilon should be chosen obtain bounds of specific tightness.
adaptive_eps	logical: if TRUE, starts with a large epsilon and then decreases it gradually whenever the progress of the algorithm requires, until the value given in the argument epsilon is reached. Usually speeds up the algorithm substantially but is an unconventional approach, and there is no guarantee that the algorithm converges appropriately towards bounds with the tightness given by the argument epsilon.
ncores	the number of cores to be used in parallel computing.
print_progress	logical: should the progress of the algorithm be printed?

Details

The upper and lower bounds are calculated using the Gripenberg's (1996) branch-and-bound method, which is also discussed in Chang and Blondel (2013). This function can be generally used for approximating the JSR of a set of square matrices, but the main intention is STVAR applications (for models created with sstvars, the function bound_JSR should be preferred). Specifically, Kheifets and Saikkonen (2020) show that if the joint spectral radius of the companion form AR matrices of the regimes is smaller than one, the STVAR process is ergodic stationary. Virolainen (2024) shows the same result for his parametrization of of threshold and smooth transition vector autoregressive models. Therefore, if the upper bound is smaller than one, the process is stationary ergodic. However, as the condition is not necessary but sufficient and also because the bound might be too conservative, upper bound larger than one does not imply that the process is not ergodic stationary. You can try higher accuracy, and if the bound is still larger than one, the result does not tell whether the process is ergodic stationary or not.

Note that with high precision (small epsilon), the computational effort required are substantial and the estimation may take long, even though the function takes use of parallel computing. This is because with small epsilon the the number of candidate solutions in each iteration may grow exponentially and a large number of iterations may be required. For this reason, adaptive_eps=TRUE can be considered for large matrices, in which case the algorithm starts with a large epsilon, and then decreases it when new candidate solutions are not found, until the epsilon given by the argument epsilon is reached.

Value

Returns an upper bound for the joint spectral radius of the "companion form AR matrices" of the regimes.

References

• C-T Chang and V.D. Blondel. 2013 . An experimental study of approximation algorithms for the joint spectral radius. Numerical algorithms, 64, 181-202.

• Gripenberg, G. 1996. Computing the joint spectral radius. Linear Algebra and its Applications, 234, 43–60.

• I.L. Kheifets, P.J. Saikkonen. 2020. Stationarity and ergodicity of Vector STAR models. Econometric Reviews, 39:4, 407-414.

• Virolainen S. 2024. Identification by non-Gaussianity in structural threshold and smooth transition vector autoregressive models. Unpublished working paper, available in ArXiv.

See Also

bound_JSR

Examples

## Below examples take approximately 5 seconds to run.

# A set of two (5x5) square matrices:
set.seed(1); S1 <- array(rnorm(20*20*2), dim=c(5, 5, 2))

# Bound the joint spectral radius of the set of matrices S1, with the
# approximate tightness epsilon=0.01:
bound_jsr_G(S1, epsilon=0.01, adaptive_eps=FALSE)

# Obtain bounds faster with adaptive_eps=TRUE:
bound_jsr_G(S1, epsilon=0.01, adaptive_eps=TRUE)
# Note that the upper bound is not the same as with adaptive_eps=FALSE.

# A set of three (3x3) square matrices:
set.seed(2); S2 <- array(rnorm(3*3*3), dim=c(3, 3, 3))

# Bound the joint spectral radius of the set of matrices S2:
bound_jsr_G(S2, epsilon=0.01, adaptive_eps=FALSE)

# Larger epsilon terminates the iteration earlier and results in wider bounds:
bound_jsr_G(S2, epsilon=0.05, adaptive_eps=FALSE)

# A set of eight (2x2) square matrices:
set.seed(3); S3 <- array(rnorm(2*2*8), dim=c(2, 2, 8))

# Bound the joint spectral radius of the set of matrices S3:
bound_jsr_G(S3, epsilon=0.01, adaptive_eps=FALSE)

---

Calculate gradient or Hessian matrix

Description

calc_gradient or calc_hessian calculates the gradient or Hessian matrix of the given function at the given point using central difference numerical approximation. get_gradient or get_hessian calculates the gradient or Hessian matrix of the log-likelihood function at the parameter estimates of a class 'stvar' object. get_soc returns eigenvalues of the Hessian matrix, and get_foc is the same as get_gradient but named conveniently.

Usage

calc_gradient(x, fn, h = 0.001, ...)

calc_hessian(x, fn, h = 0.001, ...)

get_gradient(stvar, ...)

get_hessian(stvar, ...)

get_foc(stvar, ...)

get_soc(stvar, ...)

Arguments

x	a numeric vector specifying the point where the gradient or Hessian should be calculated.
fn	a function that takes in argument x as the first argument.
h	difference used to approximate the derivatives: either a positive real number of a vector of positive real numbers with the same length as x.
...	other arguments passed to fn (or argument passed to calc_gradient or calc_hessian).
stvar	object of class "stvar"

Details

In particular, the functions get_foc and get_soc can be used to check whether the found estimates denote a (local) maximum point, a saddle point, or something else. Note that profile log-likelihood functions can be conveniently plotted with the function profile_logliks.

Value

Gradient functions return numerical approximation of the gradient and Hessian functions return numerical approximation of the Hessian. get_soc returns eigenvalues of the Hessian matrix.

Warning

No argument checks!

Examples

# Create a simple function:
foo <- function(x) x^2 + x

# Calculate the gradient at x=1 and x=-0.5:
calc_gradient(x=1, fn=foo)
calc_gradient(x=-0.5, fn=foo)

# Create a more complicated function
foo <- function(x, a, b) a*x[1]^2 - b*x[2]^2

# Calculate the gradient at x=c(1, 2) with parameter values a=0.3 and b=0.1:
calc_gradient(x=c(1, 2), fn=foo, a=0.3, b=0.1)

# Create a linear Gaussian VAR p=1 model:
theta_112 <- c(0.649526, 0.066507, 0.288526, 0.021767, -0.144024, 0.897103,
 0.601786, -0.002945, 0.067224)
mod112 <- STVAR(data=gdpdef, p=1, M=1, params=theta_112)

# Calculate the gradient of the log-likelihood function about the parameter values:
get_foc(mod112)

# Calculate the eigenvalues of the Hessian matrix of the log-likelihood function
# about the parameter values:
get_soc(mod112)

---

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

Description

check_params checks whether the parameter vector is in the parameter space.

Usage

check_params(
  data,
  p,
  M,
  d,
  params,
  weight_function = c("relative_dens", "logistic", "mlogit", "exponential", "threshold",
    "exogenous"),
  weightfun_pars = NULL,
  cond_dist = c("Gaussian", "Student", "ind_Student", "ind_skewed_t"),
  parametrization = c("intercept", "mean"),
  identification = c("reduced_form", "recursive", "heteroskedasticity",
    "non-Gaussianity"),
  AR_constraints = NULL,
  mean_constraints = NULL,
  weight_constraints = NULL,
  B_constraints = NULL,
  transition_weights,
  allow_unstab = FALSE,
  stab_tol = 0.001,
  posdef_tol = 1e-08,
  distpar_tol = 1e-08,
  weightpar_tol = 1e-08
)

Arguments

data	a matrix or class 'ts' object with d>1 columns. Each column is taken to represent a univariate time series. Missing values are not supported.
p	a positive integer specifying the autoregressive order
M	a positive integer specifying the number of regimes
d	the number of time series in the system, i.e., the dimension
params	a real valued vector specifying the parameter values. Should have the form $\theta = (\phi_{1,0},...,\phi_{M,0},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu)$, where (see exceptions below): $\phi_{m,0} =$ the $(d \times 1)$ intercept (or mean) vector of the $m$th regime. $\varphi_m = (vec(A_{m,1}),...,vec(A_{m,p}))$ $(pd^2 \times 1)$. if cond_dist="Gaussian" or "Student": $\sigma = (vech(\Omega_1),...,vech(\Omega_M))$ $(Md(d + 1)/2 \times 1)$. if cond_dist="ind_Student" or "ind_skewed_t": $\sigma = (vec(B_1),...,vec(B_M)$ $(Md^2 \times 1)$. $\alpha =$ the $(a\times 1)$ vector containing the transition weight parameters (see below). if cond_dist = "Gaussian"): Omit $\nu$ from the parameter vector. if cond_dist="Student": $\nu > 2$ is the single degrees of freedom parameter. if cond_dist="ind_Student": $\nu = (\nu_1,...,\nu_d)$ $(d \times 1)$, $\nu_i > 2$. if cond_dist="ind_skewed_t": $\nu = (\nu_1,...,\nu_d,\lambda_1,...,\lambda_d)$ $(2d \times 1)$, $\nu_i > 2$ and $\lambda_i \in (0, 1)$. For models with... weight_function="relative_dens": $\alpha = (\alpha_1,...,\alpha_{M-1})$ $(M - 1 \times 1)$, where $\alpha_m$ $(1\times 1), m=1,...,M-1$ are the transition weight parameters. weight_function="logistic": $\alpha = (c,\gamma)$ $(2 \times 1)$, where $c\in\mathbb{R}$ is the location parameter and $\gamma >0$ is the scale parameter. weight_function="mlogit": $\alpha = (\gamma_1,...,\gamma_M)$ $((M-1)k\times 1)$, where $\gamma_m$ $(k\times 1)$, $m=1,...,M-1$ contains the multinomial logit-regression coefficients of the $m$th regime. Specifically, for switching variables with indices in $I\subset\lbrace 1,...,d\rbrace$, and with $\tilde{p}\in\lbrace 1,...,p\rbrace$ lags included, $\gamma_m$ contains the coefficients for the vector $z_{t-1} = (1,\tilde{z}_{\min\lbrace I\rbrace},...,\tilde{z}_{\max\lbrace I\rbrace})$, where $\tilde{z}_{i} =(y_{it-1},...,y_{it-\tilde{p}})$, $i\in I$. So $k=1+|I|\tilde{p}$ where $|I|$ denotes the number of elements in $I$. weight_function="exponential": $\alpha = (c,\gamma)$ $(2 \times 1)$, where $c\in\mathbb{R}$ is the location parameter and $\gamma >0$ is the scale parameter. weight_function="threshold": $\alpha = (r_1,...,r_{M-1})$ $(M-1 \times 1)$, where $r_1,...,r_{M-1}$ are the thresholds. weight_function="exogenous": Omit $\alpha$ from the parameter vector. AR_constraints: Replace $\varphi_1,...,\varphi_M$ with $\psi$ as described in the argument AR_constraints. mean_constraints: Replace $\phi_{1,0},...,\phi_{M,0}$ with $(\mu_{1},...,\mu_{g})$ where $\mu_i, \ (d\times 1)$ is the mean parameter for group $i$ and $g$ is the number of groups. weight_constraints: If linear constraints are imposed, replace $\alpha$ with $\xi$ as described in the argument weigh_constraints. If weight functions parameters are imposed to be fixed values, simply drop $\alpha$ from the parameter vector. identification="heteroskedasticity": $\sigma = (vec(W),\lambda_2,...,\lambda_M)$, where $W$ $(d\times d)$ and $\lambda_m$ $(d\times 1)$, $m=2,...,M$, satisfy $\Omega_1=WW'$ and $\Omega_m=W\Lambda_mW'$, $\Lambda_m=diag(\lambda_{m1},...,\lambda_{md})$, $\lambda_{mi}>0$, $m=2,...,M$, $i=1,...,d$. B_constraints: For models identified by heteroskedasticity, replace $vec(W)$ with $\tilde{vec}(W)$ that stacks the columns of the matrix $W$ in to vector so that the elements that are constrained to zero are not included. For models identified by non-Gaussianity, replace $vec(B_1),...,vec(B_M)$ with similarly with vectorized versions $B_m$ so that the elements that are constrained to zero are not included. Above, $\phi_{m,0}$ is the intercept parameter, $A_{m,i}$ denotes the $i$th coefficient matrix of the $m$th regime, $\Omega_{m}$ denotes the positive definite error term covariance matrix of the $m$th regime, and $B_m$ is the invertible $(d\times d)$ impact matrix of the $m$th regime. $\nu_m$ is the degrees of freedom parameter of the $m$th regime. If parametrization=="mean", just replace each $\phi_{m,0}$ with regimewise mean $\mu_{m}$. $vec()$ is vectorization operator that stacks columns of a given matrix into a vector. $vech()$ stacks columns of a given matrix from the principal diagonal downwards (including elements on the diagonal) into a vector. $Bvec()$ is a vectorization operator that stacks the columns of a given impact matrix $B_m$ into a vector so that the elements that are constrained to zero by the argument B_constraints are excluded.
weight_function	What type of transition weights $\alpha_{m,t}$ should be used? "relative_dens": $\alpha_{m,t}= \frac{\alpha_mf_{m,dp}(y_{t-1},...,y_{t-p+1})}{\sum_{n=1}^M\alpha_nf_{n,dp}(y_{t-1},...,y_{t-p+1})}$, where $\alpha_m\in (0,1)$ are weight parameters that satisfy $\sum_{m=1}^M\alpha_m=1$ and $f_{m,dp}(\cdot)$ is the $dp$-dimensional stationary density of the $m$th regime corresponding to $p$ consecutive observations. Available for Gaussian conditional distribution only. "logistic": $M=2$, $\alpha_{1,t}=1-\alpha_{2,t}$, and $\alpha_{2,t}=[1+\exp\lbrace -\gamma(y_{it-j}-c) \rbrace]^{-1}$, where $y_{it-j}$ is the lag $j$ observation of the $i$th variable, $c$ is a location parameter, and $\gamma > 0$ is a scale parameter. "mlogit": $\alpha_{m,t}=\frac{\exp\lbrace \gamma_m'z_{t-1} \rbrace} {\sum_{n=1}^M\exp\lbrace \gamma_n'z_{t-1} \rbrace}$, where $\gamma_m$ are coefficient vectors, $\gamma_M=0$, and $z_{t-1}$ $(k\times 1)$ is the vector containing a constant and the (lagged) switching variables. "exponential": $M=2$, $\alpha_{1,t}=1-\alpha_{2,t}$, and $\alpha_{2,t}=1-\exp\lbrace -\gamma(y_{it-j}-c) \rbrace$, where $y_{it-j}$ is the lag $j$ observation of the $i$th variable, $c$ is a location parameter, and $\gamma > 0$ is a scale parameter. "threshold": $\alpha_{m,t} = 1$ if $r_{m-1}<y_{it-j}\leq r_{m}$ and $0$ otherwise, where $-\infty\equiv r_0<r_1<\cdots <r_{M-1}<r_M\equiv\infty$ are thresholds $y_{it-j}$ is the lag $j$ observation of the $i$th variable. "exogenous": Exogenous nonrandom transition weights, specify the weight series in weightfun_pars. See the vignette for more details about the weight functions.
weightfun_pars	If weight_function == "relative_dens": Not used. If weight_function %in% c("logistic", "exponential", "threshold"): a numeric vector with the switching variable $i\in\lbrace 1,...,d \rbrace$ in the first and the lag $j\in\lbrace 1,...,p \rbrace$ in the second element. If weight_function == "mlogit": a list of two elements: The first element $vars: a numeric vector containing the variables that should used as switching variables in the weight function in an increasing order, i.e., a vector with unique elements in $\lbrace 1,...,d \rbrace$. The second element $lags: an integer in $\lbrace 1,...,p \rbrace$ specifying the number of lags to be used in the weight function. If weight_function == "exogenous": a size (nrow(data) - p x M) matrix containing the exogenous transition weights as [t, m] for time $t$ and regime $m$. Each row needs to sum to one and only weakly positive values are allowed.
cond_dist	specifies the conditional distribution of the model as "Gaussian", "Student", "ind_Student", or "ind_skewed_t", where "ind_Student" the Student's $t$ distribution with independent components, and "ind_skewed_t" is the skewed $t$ distribution with independent components (see Hansen, 1994).
parametrization	"intercept" or "mean" determining whether the model is parametrized with intercept parameters $\phi_{m,0}$ or regime means $\mu_{m}$, m=1,...,M.
identification	is it reduced form model or an identified structural model; if the latter, how is it identified (see the vignette or the references for details)? "reduced_form": Reduced form model. "recursive": The usual lower-triangular recursive identification of the shocks via their impact responses. "heteroskedasticity": Identification by conditional heteroskedasticity, which imposes constant relative impact responses for each shock. "non-Gaussianity": Identification by non-Gaussianity; requires mutually independent non-Gaussian shocks, thus, currently available only with the conditional distribution "ind_Student".
AR_constraints	a size $(Mpd^2 \times q)$ constraint matrix $C$ specifying linear constraints to the autoregressive parameters. The constraints are of the form $(\varphi_{1},...,\varphi_{M}) = C\psi$, where $\varphi_{m} = (vec(A_{m,1}),...,vec(A_{m,p})) \ (pd^2 \times 1),\ m=1,...,M$, contains the coefficient matrices and $\psi$ $(q \times 1)$ contains the related parameters. For example, to restrict the AR-parameters to be the identical across the regimes, set $C =$ [I:...:I]' $(Mpd^2 \times pd^2)$ where I = diag(p*d^2).
mean_constraints	Restrict the mean parameters of some regimes to be identical? Provide a list of numeric vectors such that each numeric vector contains the regimes that should share the common mean parameters. For instance, if M=3, the argument list(1, 2:3) restricts the mean parameters of the second and third regime to be identical but the first regime has freely estimated (unconditional) mean. Ignore or set to NULL if mean parameters should not be restricted to be the same among any regimes. This constraint is available only for mean parametrized models; that is, when parametrization="mean".
weight_constraints	a list of two elements, $R$ in the first element and $r$ in the second element, specifying linear constraints on the transition weight parameters $\alpha$. The constraints are of the form $\alpha = R\xi + r$, where $R$ is a known $(a\times l)$ constraint matrix of full column rank ($a$ is the dimension of $\alpha$), $r$ is a known $(a\times 1)$ constant, and $\xi$ is an unknown $(l\times 1)$ parameter. Alternatively, set $R=0$ to constrain the weight parameters to the constant $r$ (in this case, $\alpha$ is dropped from the constrained parameter vector).
B_constraints	a $(d \times d)$ matrix with its entries imposing constraints on the impact matrix $B_t$: NA indicating that the element is unconstrained, a positive value indicating strict positive sign constraint, a negative value indicating strict negative sign constraint, and zero indicating that the element is constrained to zero. Currently only available for models with identification="heteroskedasticity" or "non-Gaussianity" due to the (in)availability of appropriate parametrizations that allow such constraints to be imposed.
transition_weights	(optional; only for models with cond_dist="ind_Student" or identification="non-Gaussianity") A $T \times M$ matrix containing the transition weights. If cond_dist="ind_Student" checks that the impact matrix $\sum_{m=1}^M\alpha_{m,t}^{1/2}B_m$ is invertible for all $t=1,...,T$.
allow_unstab	If TRUE, estimates not satisfying the stability condition are allowed. Always FALSE if weight_function="relative_dens".
stab_tol	numerical tolerance for stability of condition of the regimes: if the "bold A" matrix of any regime has eigenvalues larger that 1 - stat_tol the parameter is considered to be outside the parameter space. Note that if tolerance is too small, numerical evaluation of the log-likelihood might fail and cause error.
posdef_tol	numerical tolerance for positive definiteness of the error term covariance matrices: if the error term covariance matrix of any regime has eigenvalues smaller than this, the parameter is considered to be outside the parameter space. Note that if the tolerance is too small, numerical evaluation of the log-likelihood might fail and cause error.
distpar_tol	the parameter vector is considered to be outside the parameter space if the degrees of freedom parameters is not larger than 2 + distpar_tol (applies only if cond_dist="Student").
weightpar_tol	numerical tolerance for weight parameters being in the parameter space. Values closer to to the border of the parameter space than this are considered to be "outside" the parameter space.

Value

Throws an informative error if there is something wrong with the parameter vector.

References

• Kheifets I.L., Saikkonen P.J. 2020. Stationarity and ergodicity of Vector STAR models. Econometric Reviews, 39:4, 407-414.

• Lütkepohl H. 2005. New Introduction to Multiple Time Series Analysis, Springer.

• 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 models. Unpublished working paper, available as arXiv:2404.19707.

@keywords internal

Examples

# There examples will cause an informative error
 params112_notpd <- c(6.5e-01, 7.0e-01, 2.9e-01, 2.0e-02, -1.4e-01,
  9.0e-01, 6.0e-01, -1.0e-02, 1.0e-07)
 try(check_params(p=1, M=1, d=2, params=params112_notpd))

 params112_notstat <- c(6.5e-01, 7.0e-01, 10.9e-01, 2.0e-02, -1.4e-01,
  9.0e-01, 6.0e-01, -1.0e-02, 1.0e-07)
 try(check_params(p=1, M=1, d=2, params=params112_notstat))

 params112_wronglength <- c(6.5e-01, 7.0e-01, 2.9e-01, 2.0e-02, -1.4e-01,
  9.0e-01, 6.0e-01, -1.0e-02)
 try(check_params(p=1, M=1, d=2, params=params112_wronglength))

---

Simultaneously diagonalize two covariance matrices

Description

diag_Omegas Simultaneously diagonalizes two covariance matrices using eigenvalue decomposition.

Usage

diag_Omegas(Omega1, Omega2)

Arguments

Omega1	a positive definite $(dxd)$ covariance matrix $(d>1)$
Omega2	another positive definite $(dxd)$ covariance matrix

Details

See the return value and Muirhead (1982), Theorem A9.9 for details.

Value

Returns a length $d^2 + d$ vector where the first $d^2$ elements are $vec(W)$ with the columns of $W$ being (specific) eigenvectors of the matrix $\Omega_2\Omega_1^{-1}$ and the rest $d$ elements are the corresponding eigenvalues "lambdas". The result satisfies $WW' = Omega1$ and $Wdiag(lambdas)W' = Omega2$.

If Omega2 is not supplied, returns a vectorized symmetric (and pos. def.) square root matrix of Omega1.

Warning

No argument checks! Does not work with dimension $d=1$!

References

• Muirhead R.J. 1982. Aspects of Multivariate Statistical Theory, Wiley.

Examples

# Create two (2x2) coviance matrices using the parameters W and lambdas:
d <- 2 # The dimension
W0 <- matrix(1:(d^2), nrow=2) # W
lambdas0 <- 1:d # The eigenvalues
(Omg1 <- W0%*%t(W0)) # The first covariance matrix
(Omg2 <- W0%*%diag(lambdas0)%*%t(W0)) # The second covariance matrix

# Then simultaneously diagonalize the covariance matrices:
res <- diag_Omegas(Omg1, Omg2)

# Recover W:
W <- matrix(res[1:(d^2)], nrow=d, byrow=FALSE)
tcrossprod(W) # == Omg1, the first covariance matrix

# Recover lambdas:
lambdas <- res[(d^2 + 1):(d^2 + d)]
W%*%diag(lambdas)%*%t(W) # == Omg2, the second covariance matrix

---

Residual diagnostic plot for a STVAR model

Description

diagnostic_plot plots a multivariate residual diagnostic plot for either autocorrelation, conditional heteroskedasticity, or distribution, or simply draws the residual time series.

Usage

diagnostic_plot(
  stvar,
  type = c("all", "series", "ac", "ch", "dist"),
  resid_type = c("standardized", "raw"),
  maxlag = 12
)

Arguments

stvar	object of class "stvar"
type	which type of diagnostic plot should be plotted? "all" all below sequentially. "series" the residual time series. "ac" the residual autocorrelation and cross-correlation functions. "ch" the squared residual autocorrelation and cross-correlation functions. "dist" the residual histogram with theoretical density (dashed line) and QQ-plots.
resid_type	should standardized or raw residuals be used?
maxlag	the maximum lag considered in types "ac" and "ch".

Details

Auto- and cross-correlations (types "ac" and "ch") are calculated with the function acf from the package stats and the plot method for class 'acf' objects is employed.

If cond_dist == "Student" or "ind_Student", the estimates of the degrees of freedom parameters is used in theoretical densities and quantiles. If cond_dist == "ind_skewed_t", the estimates of the degrees of freedom and skewness parameters are used in theoretical densities and quantiles, and the quantile function is computed numerically.

Value

No return value, called for its side effect of plotting the diagnostic plot.

References

• Anderson H., Vahid F. 1998. Testing multiple equation systems for common nonlinear components. Journal of Econometrics, 84:1, 1-36.

• Hansen B.E. 1994. Autoregressive Conditional Density estimation. Journal of Econometrics, 35:3, 705-730.

• Kheifets I.L., Saikkonen P.J. 2020. Stationarity and ergodicity of Vector STAR models. International Economic Review, 35:3, 407-414.

• 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.

• Lütkepohl H. 2005. New Introduction to Multiple Time Series Analysis, Springer.

• McElroy T. 2017. Computation of vector ARMA autocovariances. Statistics and Probability Letters, 124, 92-96.

• Kilian L., Lütkepohl H. 20017. Structural Vector Autoregressive Analysis. 1st edition. Cambridge University Press, Cambridge.

• Tsay R. 1998. Testing and Modeling Multivariate Threshold Models. Journal of the American Statistical Association, 93:443, 1188-1202.

• Virolainen S. 2024. Identification by non-Gaussianity in structural threshold and smooth transition vector autoregressive models. Unpublished working paper, available as arXiv:2404.19707.

See Also

Portmanteau_test, profile_logliks, fitSTVAR, STVAR, LR_test, Wald_test, Rao_test

Examples

## Gaussian STVAR p=1, M=2 model, with weighted relative stationary densities
# of the regimes as the transition weight function:
theta_122relg <- c(0.734054, 0.225598, 0.705744, 0.187897, 0.259626, -0.000863,
  -0.3124, 0.505251, 0.298483, 0.030096, -0.176925, 0.838898, 0.310863, 0.007512,
  0.018244, 0.949533, -0.016941, 0.121403, 0.573269)
mod122 <- STVAR(data=gdpdef, p=1, M=2, params=theta_122relg)

# Autocorelation function of raw residuals for checking remaining autocorrelation:
diagnostic_plot(mod122, type="ac", resid_type="raw")

# Autocorelation function of squared standardized residuals for checking remaining
# conditional heteroskedasticity:
diagnostic_plot(mod122, type="ch", resid_type="standardized")

# Below, ACF of squared raw residuals, which is not very informative for evaluating
# adequacy to capture conditional heteroskedasticity, since it doesn't take into account
# the time-varying conditional covariance matrix of the model:
diagnostic_plot(mod122, type="ch", resid_type="raw")

# Similarly, below the time series of raw residuals first, and then the
# time series of standardized residuals. The latter is more informative
# for evaluating adequacy:
diagnostic_plot(mod122, type="series", resid_type="raw")
diagnostic_plot(mod122, type="series", resid_type="standardized")

# Also similarly, histogram and Q-Q plots are more informative for standardized
# residuals when evaluating model adequacy:
diagnostic_plot(mod122, type="dist", resid_type="raw") # Bad fit for GDPDEF
diagnostic_plot(mod122, type="dist", resid_type="standardized") # Good fit for GDPDEF

## Linear Gaussian VAR p=1 model:
theta_112 <- c(0.649526, 0.066507, 0.288526, 0.021767, -0.144024, 0.897103,
  0.601786, -0.002945, 0.067224)
mod112 <- STVAR(data=gdpdef, p=1, M=1, params=theta_112)
diagnostic_plot(mod112, resid_type="standardized") # All plots for std. resids
diagnostic_plot(mod112, resid_type="raw") # All plots for raw residuals

---

Filter inappropriate the estimates produced by fitSTVAR

Description

filter_estimates filters out inappropriate estimates produced by fitSTVAR: can be used to obtain the (possibly) appropriate estimate with the largest found log-likelihood (among possibly appropriate estimates) as well as (possibly) appropriate estimates based on smaller log-likelihoods.

Usage

filter_estimates(
  stvar,
  which_largest = 1,
  filter_stab = TRUE,
  calc_std_errors = FALSE
)

Arguments

stvar	a class 'stvar' object defining a structural STVAR model that is identified by heteroskedasticity or non-Gaussianity, typically created with fitSSTVAR.
which_largest	an integer at least one specifying the (possibly) appropriate estimate corresponding to which largest log-likelihood should be returned. E.g., if which_largest=2, the function will return among the estimates that it does not deem inappropriate the one that has the second largest log-likelihood.
filter_stab	Should estimates close to breaking the usual stability condition be filtered out?
calc_std_errors	should approximate standard errors be calculated?

Details

The function goes through the estimates produced by fitSTVAR and checks which estimates are deemed inappropriate. That is, estimates that are not likely solutions of interest. Specifically, solutions that incorporate a near-singular error term covariance matrix (any eigenvalue less than $0.002$), any modulus of the eigenvalues of the companion form AR matrices larger than $0.9985$ (indicating the necessary condition for stationarity is close to break), or transition weights such that they are close to zero for almost all $t$ for at least one regime. Then, among the solutions are not deemed inappropriate, it returns a STVAR models based on the estimate that has the which_largest largest log-likelihood.

The function filter_estimates is kind of a version of alt_stvar that only considers estimates that are not deemed inappropriate

Value

Returns an S3 object of class 'stvar' defining a smooth transition VAR model. The returned list contains the following components (some of which may be NULL depending on the use case):

data	The input time series data.
model	A list describing the model structure.
params	The parameters of the model.
std_errors	Approximate standard errors of the parameters, if calculated.
transition_weights	The transition weights of the model.
regime_cmeans	Conditional means of the regimes, if data is provided.
total_cmeans	Total conditional means of the model, if data is provided.
total_ccovs	Total conditional covariances of the model, if data is provided.
uncond_moments	A list of unconditional moments including regime autocovariances, variances, and means.
residuals_raw	Raw residuals, if data is provided.
residuals_std	Standardized residuals, if data is provided.
structural_shocks	Recovered structural shocks, if applicable.
loglik	Log-likelihood of the model, if data is provided.
IC	The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.
all_estimates	The parameter estimates from all estimation rounds, if applicable.
all_logliks	The log-likelihood of the estimates from all estimation rounds, if applicable.
which_converged	Indicators of which estimation rounds converged, if applicable.
which_round	Indicators of which round of optimization each estimate belongs to, if applicable.
LS_estimates	The least squares estimates of the parameters in the form $(\phi_{1,0},...,\phi_{M,0},\varphi_1,...,\varphi_M,\alpha$ (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

See Also

fitSTVAR, alt_stvar

Examples

# Fit a two-regime STVAR model with logistic transition weights and Student's t errors
 fit12 <- fitSTVAR(gdpdef, p=1, M=2, weight_function="logistic", weightfun_pars=c(2, 1),
  cond_dist="Student", nrounds=2, ncores=2, seeds=c(2, 5))
 fit12

 # Filter through inappropriate estimates and obtain the second best appropriate solution:
 fit12_2 <- filter_estimates(fit12, which_largest=2)
 fit12_2 # The same model since the two estimation rounds yielded the same estimate

---

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

Description

fitSSTVAR uses a robust method and a variable metric algorithm to estimate a structural STVAR model based on preliminary estimates from a reduced form model.

Usage

fitSSTVAR(
  stvar,
  identification = c("recursive", "heteroskedasticity", "non-Gaussianity"),
  B_constraints = NULL,
  B_pm_reg = NULL,
  B_perm = NULL,
  B_signs = NULL,
  maxit = 1000,
  maxit_robust = 1000,
  robust_method = c("Nelder-Mead", "SANN", "none"),
  print_res = TRUE,
  calc_std_errors = TRUE
)

Arguments

stvar	a an object of class 'stvar', created by, e.g., fitSTVAR, specifying a reduced form or a structural model
identification	Which identification should the structural model use? (see the vignette or the references for details) "recursive": The usual lower-triangular recursive identification of the shocks via their impact responses. "heteroskedasticity": Identification by conditional heteroskedasticity, which imposes constant relative impact responses for each shock.
B_constraints	Employ further constraints on the impact matrix? A $(d \times d)$ matrix with its entries imposing constraints on the impact matrix $B_t$: NA indicating that the element is unconstrained, a positive value indicating strict positive sign constraint, a negative value indicating strict negative sign constraint, and zero indicating that the element is constrained to zero. Currently only available for models with identification="heteroskedasticity" due to the (in)availability of appropriate parametrizations that allow such constraints to be imposed.
B_pm_reg	an integer between $1$ and $M$ specifying the regime the permutations and sign changes of $B_m$ specified in the arguments B_perm and B_signs are applied to.
B_perm	a numeric vector of length $d$ specifying the permutation of the columns of the impact matrix $B_m$ of a single regime specified in the argument B_pm_reg prior to re-estimating the model. Applicable only for models with cond_dist = "ind_Student" or "ind_skewed_t".
B_signs	a numeric vector specifying the columns of the impact matrix of a single regime specified in the argument B_pm_reg that should be multiplied by -1 prior to reordering them according to B_perm (if specified). Applicable only for models with cond_dist = "ind_Student" or "ind_skewed_t".
maxit	the maximum number of iterations in the variable metric algorithm.
maxit_robust	the maximum number of iterations on the first phase robust estimation, if employed.
robust_method	Should some robust estimation method be used in the estimation before switching to the gradient based variable metric algorithm? See details.
print_res	should summaries of estimation results be printed?
calc_std_errors	should approximate standard errors be calculated?

Details

When the structural model does not impose overidentifying constraints, it is directly obtained from the reduced form model, and estimation is not required. When overidentifying constraints are imposed, the model is estimated subject to the constraints.

Using the robust estimation method before switching to the variable metric can be useful if the initial estimates are not very close to the ML estimate of the structural model, as the variable metric algorithm (usually) converges to a nearby local maximum or saddle point. However, if the initial estimates are far from the ML estimate, the resulting solution is likely local only due to the complexity of the model. Note that Nelder-Mead algorithm is much faster than SANN but can get stuck at a local solution. This is particularly the case when the imposed overidentifying restrictions are such that the unrestricted estimate is not close to satisfying them. Nevertheless, in most practical cases, the model is just identified and estimation is not required, and often reasonable overidentifying constraints are close to the unrestricted estimate.

Employs the estimation function optim from the package stats that implements the optimization algorithms. See ?optim for the documentation on the optimization methods.

The arguments B_pm_reg, B_perm, and B_signs can be used to explore estimates based various orderings and sign changes of the columns of the impact matrices $B_m$ of specific regimes. This can be useful in the presence of weak identification with respect to the ordering or signs of the columns $B_2,...,B_M$ (see Virolainen 2024).

Value

Returns an S3 object of class 'stvar' defining a smooth transition VAR model. The returned list contains the following components (some of which may be NULL depending on the use case):

data	The input time series data.
model	A list describing the model structure.
params	The parameters of the model.
std_errors	Approximate standard errors of the parameters, if calculated.
transition_weights	The transition weights of the model.
regime_cmeans	Conditional means of the regimes, if data is provided.
total_cmeans	Total conditional means of the model, if data is provided.
total_ccovs	Total conditional covariances of the model, if data is provided.
uncond_moments	A list of unconditional moments including regime autocovariances, variances, and means.
residuals_raw	Raw residuals, if data is provided.
residuals_std	Standardized residuals, if data is provided.
structural_shocks	Recovered structural shocks, if applicable.
loglik	Log-likelihood of the model, if data is provided.
IC	The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.
all_estimates	The parameter estimates from all estimation rounds, if applicable.
all_logliks	The log-likelihood of the estimates from all estimation rounds, if applicable.
which_converged	Indicators of which estimation rounds converged, if applicable.
which_round	Indicators of which round of optimization each estimate belongs to, if applicable.
LS_estimates	The least squares estimates of the parameters in the form $(\phi_{1,0},...,\phi_{M,0},\varphi_1,...,\varphi_M,\alpha$ (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

References

• Kilian L., Lütkepohl H. 20017. Structural Vector Autoregressive Analysis. 1st edition. Cambridge University Press, Cambridge.

• Lütkepohl H., Netšunajev A. 2017. Structural vector autoregressions with smooth transition in variances. Journal of Economic Dynamics & Control, 84, 43-57.

• Virolainen S. 2024. Identification by non-Gaussianity in structural threshold and smooth transition vector autoregressive models. Unpublished working paper, available as arXiv:2404.19707.

See Also

fitSTVAR, STVAR, optim

Examples

## These are long running examples that take approximately 1 minute to run.

## Estimate first a reduced form Gaussian STVAR p=3, M=2 model with the weighted relative
# stationary densities of the regimes as the transition weight function, and the means and
# AR matrices constrained to be identical across the regimes:
fit32cm <- fitSTVAR(gdpdef, p=3, M=2, AR_constraints=rbind(diag(3*2^2), diag(3*2^2)),
  weight_function="relative_dens", mean_constraints=list(1:2), parametrization="mean",
  nrounds=1, seeds=1, ncores=1)

# Then, we estimate/create various structural models based on the reduced form model.
# Create a structural model with the shocks identified recursively:
fit32cms_rec <- fitSSTVAR(fit32cm, identification="recursive")

# Create a structural model with the shocks identified by conditional heteroskedasticity:
fit32cms_hetsked <- fitSSTVAR(fit32cm, identification="heteroskedasticity")
fit32cms_hetsked # Print the estimates

# Estimate a structural model with the shocks identified by conditional heteroskedasticity
# and overidentifying constraints imposed on the impact matrix: positive diagonal element
# and zero upper right element:
fit32cms_hs2 <- fitSSTVAR(fit32cm, identification="heteroskedasticity",
 B_constraints=matrix(c(1, NA, 0, 1), nrow=2))

# Estimate a structural model with the shocks identified by conditional heteroskedasticity
# and overidentifying constraints imposed on the impact matrix: positive diagonal element
# and zero off-diagonal elements:
fit32cms_hs3 <- fitSSTVAR(fit32cms_hs2, identification="heteroskedasticity",
 B_constraints=matrix(c(1, 0, 0, 1), nrow=2))

# Estimate first a reduced form two-regime Threshold VAR p=1 model with
# with independent skewed t shocks, and the first lag of the second variable
# as the switching variable, and AR matrices constrained to be identical
# across the regimes:
fit12c <- fitSTVAR(gdpdef, p=1, M=2, cond_dist="ind_skewed_t",
 AR_constraints=rbind(diag(1*2^2), diag(1*2^2)), weight_function="threshold",
 weightfun_pars=c(2, 1), nrounds=1, seeds=1, ncores=1)

# Due to the independent non-Gaussian shocks, the structural shocks are readily
# identified. The following returns the same model but marked as structural
# with the shocks identified by non-Gaussianity:
fit12c <- fitSSTVAR(fit12c)

# Estimate a model based on a reversed ordering of the columns of the impact matrix B_2:
fit12c2 <- fitSSTVAR(fit12c, B_pm_reg=2, B_perm=c(2, 1))

# Estimate a model based on reversed signs of the second column of B_2 and reversed
# ordering of the columns of B_2:
fit12c3 <- fitSSTVAR(fit12c, B_pm_reg=2, B_perm=c(2, 1), B_signs=2)

---

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

Description

fitSTVAR estimates a reduced form smooth transition VAR model in two phases or three phases. In the two-phase procedure: in the first phase, it uses a genetic algorithm (GA) to find starting values for a gradient based variable metric algorithm (VA), which it then uses to finalize the estimation in the second phase. Parallel computing is utilized to perform multiple rounds of estimations in parallel. In the three-phase procedure, the autoregressive and weight parameters are first estimated by least squares (LS) to obtain initial estimates for GA, and the rest of the procedure proceeds as in the two-phase procedure.

Usage

fitSTVAR(
  data,
  p,
  M,
  weight_function = c("relative_dens", "logistic", "mlogit", "exponential", "threshold",
    "exogenous"),
  weightfun_pars = NULL,
  cond_dist = c("Gaussian", "Student", "ind_Student", "ind_skewed_t"),
  parametrization = c("intercept", "mean"),
  AR_constraints = NULL,
  mean_constraints = NULL,
  weight_constraints = NULL,
  estim_method,
  penalized,
  penalty_params = c(0.05, 0.2),
  allow_unstab,
  nrounds,
  ncores = 2,
  maxit = 1000,
  seeds = NULL,
  print_res = TRUE,
  use_parallel = TRUE,
  calc_std_errors = TRUE,
  ...
)

Arguments

data	a matrix or class 'ts' object with d>1 columns. Each column is taken to represent a univariate time series. Missing values are not supported.
p	a positive integer specifying the autoregressive order
M	a positive integer specifying the number of regimes
weight_function	What type of transition weights $\alpha_{m,t}$ should be used? "relative_dens": $\alpha_{m,t}= \frac{\alpha_mf_{m,dp}(y_{t-1},...,y_{t-p+1})}{\sum_{n=1}^M\alpha_nf_{n,dp}(y_{t-1},...,y_{t-p+1})}$, where $\alpha_m\in (0,1)$ are weight parameters that satisfy $\sum_{m=1}^M\alpha_m=1$ and $f_{m,dp}(\cdot)$ is the $dp$-dimensional stationary density of the $m$th regime corresponding to $p$ consecutive observations. Available for Gaussian conditional distribution only. "logistic": $M=2$, $\alpha_{1,t}=1-\alpha_{2,t}$, and $\alpha_{2,t}=[1+\exp\lbrace -\gamma(y_{it-j}-c) \rbrace]^{-1}$, where $y_{it-j}$ is the lag $j$ observation of the $i$th variable, $c$ is a location parameter, and $\gamma > 0$ is a scale parameter. "mlogit": $\alpha_{m,t}=\frac{\exp\lbrace \gamma_m'z_{t-1} \rbrace} {\sum_{n=1}^M\exp\lbrace \gamma_n'z_{t-1} \rbrace}$, where $\gamma_m$ are coefficient vectors, $\gamma_M=0$, and $z_{t-1}$ $(k\times 1)$ is the vector containing a constant and the (lagged) switching variables. "exponential": $M=2$, $\alpha_{1,t}=1-\alpha_{2,t}$, and $\alpha_{2,t}=1-\exp\lbrace -\gamma(y_{it-j}-c) \rbrace$, where $y_{it-j}$ is the lag $j$ observation of the $i$th variable, $c$ is a location parameter, and $\gamma > 0$ is a scale parameter. "threshold": $\alpha_{m,t} = 1$ if $r_{m-1}<y_{it-j}\leq r_{m}$ and $0$ otherwise, where $-\infty\equiv r_0<r_1<\cdots <r_{M-1}<r_M\equiv\infty$ are thresholds $y_{it-j}$ is the lag $j$ observation of the $i$th variable. "exogenous": Exogenous nonrandom transition weights, specify the weight series in weightfun_pars. See the vignette for more details about the weight functions.
weightfun_pars	If weight_function == "relative_dens": Not used. If weight_function %in% c("logistic", "exponential", "threshold"): a numeric vector with the switching variable $i\in\lbrace 1,...,d \rbrace$ in the first and the lag $j\in\lbrace 1,...,p \rbrace$ in the second element. If weight_function == "mlogit": a list of two elements: The first element $vars: a numeric vector containing the variables that should used as switching variables in the weight function in an increasing order, i.e., a vector with unique elements in $\lbrace 1,...,d \rbrace$. The second element $lags: an integer in $\lbrace 1,...,p \rbrace$ specifying the number of lags to be used in the weight function. If weight_function == "exogenous": a size (nrow(data) - p x M) matrix containing the exogenous transition weights as [t, m] for time $t$ and regime $m$. Each row needs to sum to one and only weakly positive values are allowed.
cond_dist	specifies the conditional distribution of the model as "Gaussian", "Student", "ind_Student", or "ind_skewed_t", where "ind_Student" the Student's $t$ distribution with independent components, and "ind_skewed_t" is the skewed $t$ distribution with independent components (see Hansen, 1994).
parametrization	"intercept" or "mean" determining whether the model is parametrized with intercept parameters $\phi_{m,0}$ or regime means $\mu_{m}$, m=1,...,M.
AR_constraints	a size $(Mpd^2 \times q)$ constraint matrix $C$ specifying linear constraints to the autoregressive parameters. The constraints are of the form $(\varphi_{1},...,\varphi_{M}) = C\psi$, where $\varphi_{m} = (vec(A_{m,1}),...,vec(A_{m,p})) \ (pd^2 \times 1),\ m=1,...,M$, contains the coefficient matrices and $\psi$ $(q \times 1)$ contains the related parameters. For example, to restrict the AR-parameters to be the identical across the regimes, set $C =$ [I:...:I]' $(Mpd^2 \times pd^2)$ where I = diag(p*d^2).
mean_constraints	Restrict the mean parameters of some regimes to be identical? Provide a list of numeric vectors such that each numeric vector contains the regimes that should share the common mean parameters. For instance, if M=3, the argument list(1, 2:3) restricts the mean parameters of the second and third regime to be identical but the first regime has freely estimated (unconditional) mean. Ignore or set to NULL if mean parameters should not be restricted to be the same among any regimes. This constraint is available only for mean parametrized models; that is, when parametrization="mean".
weight_constraints	a list of two elements, $R$ in the first element and $r$ in the second element, specifying linear constraints on the transition weight parameters $\alpha$. The constraints are of the form $\alpha = R\xi + r$, where $R$ is a known $(a\times l)$ constraint matrix of full column rank ($a$ is the dimension of $\alpha$), $r$ is a known $(a\times 1)$ constant, and $\xi$ is an unknown $(l\times 1)$ parameter. Alternatively, set $R=0$ to constrain the weight parameters to the constant $r$ (in this case, $\alpha$ is dropped from the constrained parameter vector).
estim_method	either "two-phase" or "three-phase" (the latter is the default option for threshold models and the former is currently the only option for other models). See details.
penalized	should penalized log-likelihood function be used that penalizes the log-likelihood function when the parameter values are close the boundary of the stability region or outside it? If TRUE, estimates that do not satisfy the stability condition are allowed (except when weight_function="relative_dens"). The default is TRUE for three-phase estimation and FALSE for two-phase estimation.
penalty_params	a numeric vector with two positive elements specifying the penalization parameters: the first element determined how far from the boundary of the stability region the penalization starts (a number between zero and one, smaller number starts penalization closer to the boundary) and the second element is a tuning parameter for the penalization (a positive real number, a higher value penalizes non-stability more).
allow_unstab	If TRUE, estimates not satisfying the stability condition are allowed. Always FALSE if weight_function="relative_dens".
nrounds	the number of estimation rounds that should be performed. The default is (M*ncol(data))^3 when estim_method="two-phase" and (M*ncol(data))^2 when estim_method="three-phase".
ncores	the number CPU cores to be used in parallel computing.
maxit	the maximum number of iterations in the variable metric algorithm.
seeds	a length nrounds vector containing the random number generator seed for each call to the genetic algorithm, or NULL for not initializing the seed.
print_res	should summaries of estimation results be printed?
use_parallel	employ parallel computing? If use_parallel=FALSE && print_res=FALSE, nothing is printed during the estimation process.
...	additional settings passed to the function GAfit employing the genetic algorithm.

Details

If you wish to estimate a structural model, estimate first the reduced form model and then use the use the function fitSSTVAR to create (and estimate if necessary) the structural model based on the estimated reduced form model.

three-phase estimation. With estim_method="three-phase" (currently only available for threshold VAR models), an extra phase is added to the beginning of the two-phase estimation procedure: the autoregressive and weight function parameters are first estimated by the method of least squares. Then, these initial estimates are used to create an initial population to the genetic algorithm, and the rest of the procedure proceeds as in the the two-phase procedure. This allows to use substantially decrease the required number of estimation rounds, and thereby speeds up the estimation substantially. On the other hand, the three-phase procedure tends to produce estimates close to the initial LS estimates, while the two-phase procedure explores the parameter space more thoroughly.

The rest concerns both two-phase and three-phase procedures.\ Because of complexity and high multimodality of the log-likelihood function, it is not certain that the estimation algorithm will end up in the global maximum point. When estim_method="two-phase", it is expected that many of the estimation rounds will end up in some local maximum or a saddle point instead. Therefore, a (sometimes very large) number of estimation rounds is required for reliable results (when estim_method="three-phase" substantially smaller number should be sufficient). Due to identification problems and high complexity of the surface of the log-likelihood function, the estimation may fail especially in the cases where the number of regimes is chosen too large.

The estimation process is computationally heavy and it might take considerably long time for large models to estimate, particularly if estim_method="two-phase". Note that reliable estimation of model with cond_dist == "ind_Student" or "ind_skewed_t" is more difficult than with Gaussian or Student's t models due to the increased complexity.

If the iteration limit maxit in the variable metric algorithm is reached, one can continue the estimation by iterating more with the function iterate_more. Alternatively, one may use the found estimates as starting values for the genetic algorithm and employ another round of estimation (see ??GAfit how to set up an initial population with the dot parameters).

If the estimation algorithm performs poorly, it usually helps to scale the individual series so that they vary roughly in the same scale. This makes it is easier to draw reasonable AR coefficients and (with some weight functions) weight parameter values in the genetic algorithm. Even if the estimation algorithm somewhat works, it should be preferred to scale the data so that most of the AR coefficients will not be very large, as the genetic algorithm works better with relatively small AR coefficients. If needed, another package can be used to fit linear VARs to the series to see which scaling of the series results in relatively small AR coefficients. You should avoid very small (or very high) variance in the data as well, so that the eigenvalues of the covariance matrices are in a reasonable range.

weight_constraints: If you are using weight constraints other than restricting some of the weight parameters to known constants, make sure the constraints are sensible. Otherwise, the estimation may fail due to the estimation algorithm not being able to generate reasonable random guesses for the values of the constrained weight parameters.

Filtering inappropriate estimates: fitSTVAR automatically filters through estimates that it deems "inappropriate". That is, estimates that are not likely solutions of interest. Specifically, solutions that incorporate a near-singular error term covariance matrix (any eigenvalue less than $0.002$), any modulus of the eigenvalues of the companion form AR matrices larger than $0.9985$ (indicating the necessary condition for stationarity is close to break), or transition weights such that they are close to zero for almost all $t$ for at least one regime. You can also always find the solutions of interest yourself by using the function alt_stvar as well since results from all estimation rounds are saved).

Value

Returns an S3 object of class 'stvar' defining a smooth transition VAR model. The returned list contains the following components (some of which may be NULL depending on the use case):

data	The input time series data.
model	A list describing the model structure.
params	The parameters of the model.
std_errors	Approximate standard errors of the parameters, if calculated.
transition_weights	The transition weights of the model.
regime_cmeans	Conditional means of the regimes, if data is provided.
total_cmeans	Total conditional means of the model, if data is provided.
total_ccovs	Total conditional covariances of the model, if data is provided.
uncond_moments	A list of unconditional moments including regime autocovariances, variances, and means.
residuals_raw	Raw residuals, if data is provided.
residuals_std	Standardized residuals, if data is provided.
structural_shocks	Recovered structural shocks, if applicable.
loglik	Log-likelihood of the model, if data is provided.
IC	The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.
all_estimates	The parameter estimates from all estimation rounds, if applicable.
all_logliks	The log-likelihood of the estimates from all estimation rounds, if applicable.
which_converged	Indicators of which estimation rounds converged, if applicable.
which_round	Indicators of which round of optimization each estimate belongs to, if applicable.
LS_estimates	The least squares estimates of the parameters in the form $(\phi_{1,0},...,\phi_{M,0},\varphi_1,...,\varphi_M,\alpha$ (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

S3 methods

The following S3 methods are supported for class 'stvar': logLik, residuals, print, summary, predict, simulate, and plot.

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.

• Koivisto T., Luoto J., Virolainen S. 2025. Unpublished working paper.

• 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.

• Kheifets I.L., Saikkonen P.J. 2020. Stationarity and ergodicity of Vector STAR models. Econometric Reviews, 39:4, 407-414.

• Tsay R. 1998. Testing and Modeling Multivariate Threshold Models. Journal of the American Statistical Association, 93:443, 1188-1202.

• Virolainen S. 2024. Identification by non-Gaussianity in structural threshold and smooth transition vector autoregressive models. Unpublished working paper, available as arXiv:2404.19707.

See Also

fitSSTVAR, STVAR, GAfit, iterate_more, filter_estimates

Examples

## These are long running examples. Running all the below examples will take
## approximately three minutes.
# When estimating the models in empirical applications, typically a large number
# of estimation rounds (set by the argument 'nrounds') should be used. These examples
# use only a small number of rounds to make the running time of the examples reasonable.

# The below examples make use of the two-variate dataset 'gdpdef' containing
# the the quarterly U.S. GDP and GDP deflator from 1947Q1 to 2019Q4.

# Estimate Gaussian STVAR model of autoregressive order p=3 and two regimes (M=2),
# with the weighted relative stationary densities of the regimes as the transition
# weight function. The estimation is performed with 2 rounds and 2 CPU cores, with
# the random number generator seeds set for reproducibility.
fit32 <- fitSTVAR(gdpdef, p=3, M=2, weight_function="relative_dens",
 cond_dist="Gaussian", nrounds=2, ncores=2, seeds=1:2)

# Examine the results:
fit32 # Printout of the estimates
summary(fit32) # A more detailed summary printout
plot(fit32) # Plot the fitted transition weights
get_foc(fit32) # Gradient of the log-likelihood function about the estimate
get_soc(fit32) # Eigenvalues of the Hessian of the log-lik. fn. about the estimate
profile_logliks(fit32) # Profile log-likelihood functions about the estimate

# Estimate a two-regime Student's t STVAR p=3 model with logistic transition weights
# and the first lag of the second variable as the switching variable, only two
# estimation rounds using two CPU cores:
fitlogistict32 <- fitSTVAR(gdpdef, p=3, M=2, weight_function="logistic", weightfun_pars=c(2, 1),
 cond_dist="Student", nrounds=2, ncores=2, seeds=1:2)
summary(fitlogistict32) # Summary printout of the estimates

# Estimate a two-regime threshold VAR p=3 model with independent skewed t shocks
# using the three-phase estimation procedure.
# The first lag of the the second variable is specified as the switching variable,
# and the threshold parameter constrained to the fixed value 1.
fitthres32wit <- fitSTVAR(gdpdef, p=3, M=2, weight_function="threshold", weightfun_pars=c(2, 1),
  cond_dist="ind_skewed_t", weight_constraints=list(R=0, r=1), estim_method="three-phase",
  nrounds=2, ncores=2, seeds=1:2)
plot(fitthres32wit) # Plot the fitted transition weights

# Estimate a two-regime STVAR p=1 model with exogenous transition weights defined as the indicator
# of NBER based U.S. recessions (source: St. Louis Fed database). Moreover, restrict the AR matrices
# to be identical across the regimes (i.e., allowing for time-variation in the intercepts and the
# covariance matrix only):

# Step 1: Define transition weights of Regime 1 as the indicator of NBER based U.S. recessions
# (the start date of weights is start of data + p, since the first p values are used as the initial
# values):
tw1 <- c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

# Step 2: Define the transition weights of Regime 2 as one minus the weights of Regime 1, and
# combine the weights to matrix of transition weights:
twmat <- cbind(tw1, 1 - tw1)

# Step 3: Create the appropriate constraint matrix:
C_122 <- rbind(diag(1*2^2), diag(1*2^2))

# Step 4: Estimate the model by specifying the weights in the argument 'weightfun_pars'
# and the constraint matrix in the argument 'AR_constraints':
fitexo12cit <- fitSTVAR(gdpdef, p=1, M=2, weight_function="exogenous", weightfun_pars=twmat,
  cond_dist="ind_Student", AR_constraints=C_122, nrounds=2, ncores=2, seeds=1:2)
plot(fitexo12cit) # Plot the transition weights
summary(fitexo12cit) # Summary printout of the estimates

# Estimate a two-regime Gaussian STVAR p=1 model with the weighted relative stationary densities
# of the regimes as the transition weight function, and the means of the regimes
# and AR matrices constrained to be identical across the regimes (i.e., allowing for time-varying
# conditional covariance matrix only):
fit12cm <- fitSTVAR(gdpdef, p=1, M=2, weight_function="relative_dens", cond_dist="Gaussian",
 AR_constraints=C_122, mean_constraints=list(1:2), parametrization="mean", nrounds=2, seeds=1:2)
fit12cm # Print the estimates

# Estimate a two-regime Gaussian STVAR p=1 model with the weighted relative stationary densities
# of the regimes as the transition weight function; constrain AR matrices to be identical
# across the regimes and also constrain the off-diagonal elements of the AR matrices to be zero.
mat0 <- matrix(c(1, rep(0, 10), 1, rep(0, 8), 1, rep(0, 10), 1), nrow=2*2^2, byrow=FALSE)
C_222 <- rbind(mat0, mat0) # The constraint matrix
fit22c <- fitSTVAR(gdpdef, p=2, M=2, weight_function="relative_dens", cond_dist="Gaussian",
 AR_constraints=C_222, nrounds=2, seeds=1:2)
fit22c # Print the estimates

# Estimate a two-regime Student's t STVAR p=3 model with logistic transition weights
# and the first lag of the second variable as the switching variable. Constraint the location
# parameter to the fixed value 1 and leave the scale parameter unconstrained.
fitlogistic32w <- fitSTVAR(gdpdef, p=3, M=2, weight_function="logistic", weightfun_pars=c(2, 1),
 weight_constraints=list(R=matrix(c(0, 1), nrow=2), r=c(1, 0)), nrounds=2, seeds=1:2)
plot(fitlogistic32w) # Plot the fitted transition weights

# Estimate a two-regime Gaussian STVAR p=3 model with multinomial logit transition weights
# using the second variable is the switching variable with two lags. Constrain the AR matrices
# identical across the regimes (allowing for time-variation in the intercepts and covariance
# matrix).
C_322 <- rbind(diag(3*2^2), diag(3*2^2)) # The constraint matrix
fitmlogit32c <- fitSTVAR(gdpdef, p=3, M=2, weight_function="mlogit", cond_dist="Gaussian",
 weightfun_pars=list(vars=2, lags=2), AR_constraints=C_322, nrounds=1, seeds=3, ncores=1)
plot(fitmlogit32c) # Plot the fitted transition weights

# Estimate a two-regime Gaussian STVAR p=3 model with exponential transition weights and the first
# lag of the second variable as switching variable, and AR parameter constrained identical across
# the regimes, means constrained identical across the regimes, and the location parameter
# constrained to 0.5 (but scale parameter unconstrained).
fitexp32cmw <- fitSTVAR(gdpdef, p=3, M=2, weight_function="exponential", weightfun_pars=c(2, 1),
 cond_dist="Student", AR_constraints=C_322, mean_constraints=list(1:2),
 weight_constraints=list(R=matrix(c(0, 1), nrow=2), r=c(0.5, 0)), nrounds=1, seeds=1, ncores=1)
summary(fitexp32cmw) # Summary printout of the estimates

---