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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

a positive integer specifying the number of generations to be ran through in the genetic algorithm.

a positive even integer specifying the population size in the genetic algorithm. Default is 10*n_params.

a positive integer specifying the generation after which the random mutations in the genetic algorithm are "smart". This means that mutating individuals will mostly mutate fairly close (or partially close) to the best fitting individual (which has the least regimes with time varying mixing weights practically at zero) so far.

a list of parameter vectors from which the initial population of the genetic algorithm will be generated from. The parameter vectors hould have the form \theta = (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu), where (see exceptions below):

• \phi_{m} = the (d \times 1) intercept (or mean) vector of the mth 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 mth 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 threshold values.

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},...,\phi_{M} 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.

Above, \phi_{m} is the intercept parameter, A_{m,i} denotes the ith coefficient matrix of the mth regime, \Omega_{m} denotes the positive definite error term covariance matrix of the mth regime, and B_m is the invertible (d\times d) impact matrix of the mth regime. \nu_m is the degrees of freedom parameter of the mth regime. If parametrization=="mean", just replace each \phi_{m} 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.

a size (dx1) vector defining means of the normal distributions from which each mean parameter \mu_{m} is drawn from in random mutations. Default is colMeans(data). Note that mean-parametrization is always used for optimization in GAfit - even when parametrization=="intercept". However, input (in initpop) and output (return value) parameter vectors can be intercept-parametrized.

a size (dx1) strictly positive vector defining standard deviations of the normal distributions from which each mean parameter \mu_{m} is drawn from in random mutations. Default is vapply(1:d, function(i1) sd(data[,i1]), numeric(1)).

a size (dx1) strictly positive vector specifying the scale and variability of the random covariance matrices in random mutations. The covariance matrices are drawn from (scaled) Wishart distribution. Expected values of the random covariance matrices are diag(omega_scale). Standard deviations of the diagonal elements are sqrt(2/d)*omega_scale[i] and for non-diagonal elements they are sqrt(1/d*omega_scale[i]*omega_scale[j]). Note that for d>4 this scale may need to be chosen carefully. Default in GAfit is var(stats::ar(data[,i], order.max=10)$resid, na.rm=TRUE), i=1,...,d. This argument is ignored if cond_dist == "ind_Student".

a size (d \times 1) strictly positive vector specifying the mean and variability of the random impact matrices in random mutations. In Regime 1, the mean of the error term covariance matrix implied by the random impact matrix will be 0.95*diag(B_scale) and in the rest of the regimes diag(B_scale), whereas the variability increases with B_scale. Default in GAfit is var(stats::ar(data[,i], order.max=10)$resid, na.rm=TRUE), i=1,...,d. This argument is ignored if cond_dist != "ind_Student".

For...

weight_function %in% c("relative_dens", "exogenous"):

not used.

weight_function %in% c("logistic", "exponential"):

length three vector with the mean (in the first element) and standard deviation (in the second element) of the normal distribution the location parameter is drawn from in random mutations. The third element is the standard deviation of the normal distribution from whose absolute value the location parameter is drawn from.

weight_function == "mlogit":

length two vector with the mean (in the first element) and standard deviation (in the second element) of the normal distribution the coefficients of the logit sub model's constant terms are drawn from in random mutations. The third element is the standard deviation of the normal distribution from which the non-constant regressors' coefficients are drawn from.

weight_function == "threshold":

a lenght two vector with the lower bound, in the first element and the upper bound, in the second element, of the uniform distribution threshold parameters are drawn from in random mutations.

a positive real number between zero and one adjusting how large AR parameter values are typically proposed in construction of the initial population: larger value implies larger coefficients (in absolute value). After construction of the initial population, a new scale is drawn from (0, upper_ar_scale) uniform distribution in each iteration.

the upper bound for ar_scale parameter (see above) in the random mutations. Setting this too high might lead to failure in proposing new parameters that are well enough inside the parameter space, and especially with large p one might want to try smaller upper bound (e.g., 0.5). With large p or d, upper_ar_scale is restricted from above, see the details section.

a positive real number adjusting how large AR parameter values are typically proposed in some random mutations (if AR constraints are employed, in all random mutations): larger value implies smaller coefficients (in absolute value). Values larger than 1 can be used if the AR coefficients are expected to be very small. If set smaller than 1, be careful as it might lead to failure in the creation of parameter candidates that satisfy the stability condition.

a non-negative real number specifying how much should natural selection favor individuals with less regimes that have almost all mixing weights (practically) at zero. Set to zero for no favoring or large number for heavy favoring. Without any favoring the genetic algorithm gets more often stuck in an area of the parameter space where some regimes are wasted, but with too much favouring the best genes might never mix into the population and the algorithm might converge poorly. Default is 1 and it gives 2x larger surviving probability weights for individuals with no wasted regimes compared to individuals with one wasted regime. Number 2 would give 3x larger probability weights etc.

Perform penalized LS estimation that minimizes penalized RSS in which estimates close to breaking or not satisfying the usual stability condition are penalized? If TRUE, the tuning parameter is set by the argument penalty_params[2], and the penalization starts when the eigenvalues of the companion form AR matrix are larger than 1 - penalty_params[1].

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

If TRUE, estimates not satisfying the stability condition are allowed. Always FALSE if weight_function="relative_dens".

a length 2 numeric vector specifying the criteria that is used to determine whether a regime is redundant (or "wasted") or not. Any regime m which satisfies sum(transitionWeights[,m] > red_criteria[1]) < red_criteria[2]*n_obs will be considered "redundant". One should be careful when adjusting this argument (set c(0, 0) to fully disable the 'redundant regime' features from the algorithm).

should the parameter space be constrained to areas where the transition weights do allocate enough weights to each regime compared to the number of observations in the regime? See the argument min_obs_coef_ga, and if three-phase estimation is used, carefully adjust it to be small enough so that estimates obtained from LS/NLS estimation are not discarded. A smaller number than min_obs_coef of LS/NLS estimation should be used.

if bound_by_weights=TRUE, determines the smallest accepted number of observations (times variables) from each regime relative to the number of parameters in the regime. For models with AR constraints, the number of AR matrix parameters in each regimes is simply assumed to be ncol(AR_constraints)/M. See the source code of the function loglikelihood for further details.

A number in [0,1] giving a probability of a "smart mutation" occuring randomly in each iteration before the iteration given by the argument smart_mu.

should the genetic algorithm return the best fitting individual which has "positive enough" mixing weights for as many regimes as possible ("alt_ind") or the individual which has the highest log-likelihood in general ("best_ind") but might have more wasted regimes?

a real number defining the minimum value of the log-likelihood function that will be considered. Values smaller than this will be treated as they were minval and the corresponding individuals will never survive. The default is -(10^(ceiling(log10(n_obs)) + d) - 1).

a vector containing fixed parameter values for intercept, autoregressive, and weight parameters that should be fixed in the initial population. Should have the form: (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\alpha, where

• (\phi_{m} = the (d \times 1) intercept vector of the mth regime.

• \varphi_m = (vec(A_{m,1}),...,vec(A_{m,p})) (pd^2 \times 1).

• \alpha vector of the weight parameters.

For models with...

AR_constraints:

Replace \varphi_1,...,\varphi_M with \psi as described in the argument AR_constraints.

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.

Note that fixed_params should always be in the intercept parametrization (and parametrization="intercept" should always be used). Passing this argument from fitSTVAR in does not do anything, as it is designed to be used with the three-phase estimation procedure only. Also, this argument does not do anything if the initial population is specified in the argument initpop.

should the fixed parameters be fixed in the smart mutation phase as well? If TRUE, the fixed parameters stay fixed throughout the whole GA estimation.

a single value, interpreted as an integer, or NULL, that sets seed for the random number generator in the beginning of the function call. If calling GAfit from fitSTVAR, use the argument seeds instead of passing the argument seed.

an object of class 'stvar', created by, e.g., fitSTVAR or fitSSTVAR.

a positive integer specifying the horizon how far ahead should the GFEVD be calculated.

What sign and size should be used for all shocks? By the normalization, the conditional covariance matrix of the structural error is an identity matrix.

What type initial values are used for estimating the GIRFs that the GFEVD is based on?

"data":

Estimate the GIRF for all the possible length p histories in the data.

"random":

Estimate the GIRF for several random initial values generated from the a specific regime specified by the argument init_regimes. The number of initial values is set with the argument R2.

"fixed":

Estimate the GIRF for the initial values specified with the argument init_values.

a numeric vector with values in 1,...,d (d=ncol(data)) specifying which the variables for which the impulse responses should be cumulative. Default is none.

set TRUE (recommended) for a special feature in which for every possible length p history in the data, or a subset of them if so specified in the argument data_gfevd_pars, the GFEVD is estimated for a shock that has the sign and size of the corresponding structural shock recovered from the data. See the details section.

a length two numeric vector with the following elements determining settings for initval_type="data" and use_data_shocks=TRUE:

1. An integer between 0 and M determining the (dominant) regime for which the GFEVD should be calculated (0 for all regimes).

2. A number between 0.5 and 1 determining how large transition weight a regime should have to be considered dominant in a given time period (i.e., determining which histories are used to calculate the GFEVD if the first element is not 0).

an integer in 1,...,M specifying the regime from which the initial values should be generated from (using a simulation procedure with a burn-in period). For models with Gaussian conditional distribution, it is also possible to generate the starting values from the stationary distribution of a regime. See the details section.

a size [p, d, R2] array specifying the initial values in each slice for each Monte Carlo repetition, where d is the number of time series in the system and R2 is an argument of this function. In each slice, the last row will be used as initial values for the first lag, the second last row for second lag etc. If not specified, initial values will be drawn from the regime specified in init_regimes.

the number of repetitions used to estimate GIRF for each initial value.

the number of initial values to be drawn/used if initval_type="random" or "fixed".

the number CPU cores to be used in parallel computing. Only single core computing is supported if an initial value is specified (and the GIRF won't thus be estimated multiple times).

Burn-in period for simulating initial values from a regime.

if weight_function="exogenous", provide a size (N+1 \times M) matrix of exogenous transition weights for the regimes: [h, m] for the (after-the-impact) period h-1 and regime m weight ([1, m] is for the impact period). Ignored if weight_function != "exogenous".

a numeric vector containing the random number generator seed for estimation of each GIRF. Should have the length of at least...

• ...nrow(data) - p + 1 if initval_type="data" or use_data_shocks=TRUE.

• ...R2 if initval_type="random".

• ...1 if initval_type="fixed.".

If the length of the vector is greater than what is needed, the extra seeds are dropped from the end of the seeds vector. Set to NULL for not initializing the seed.

employ parallel computing? If FALSE, does not print anything.

object of class 'gfevd' generated by the function GFEVD.

graphical parameters passed to the 'ts' plot method when using data_shock_pars.

if use_data_shocks, alternative plot method can be used that plots the relative contribution of a given shock to the forecast error variance of each variable at a given horizon. Should be a length two numeric vector with the number of the shock (1,..,d) in the first element and the horizon (0,1,2,...,N) in the second element.

the number of decimals to print

an integer specifying the horizon how far to print the estimates. The default is that all the values are printed.

an object of class 'stvar', created by, e.g., fitSTVAR or fitSSTVAR.

a numeric vector of length at most d (=ncol(data)) and elements in 1,...,d specifying the structural shocks for which the GIRF should be estimated.

a non-zero scalar value specifying the common size for all scalar components of the structural shock. Note that the conditional covariance matrix of the structural shock is normalized to an identity matrix and that the (generalized) impulse responses may not be symmetric with respect to the sign and size of the shock.

a positive integer specifying the horizon how far ahead should the generalized impulse responses be calculated.

the number of repetitions used to estimate GIRF for each initial value.

the number of initial values to use, i.e., to draw from init_regime if init_values are not specified. The confidence bounds will be sample quantiles of the GIRFs based on different initial values. Ignored if the argument init_value is specified. @param init_regime an integer in 1,...,M specifying the regime from which the initial values should be generated from (see ?simulate.stvar). If use_data_shocks=TRUE this is argument not used and data_girf_pars should be specified instead.

an integer in 1,...,M specifying the regime from which the initial values should be generated from (using a simulation procedure with a burn-in period). For models with Gaussian conditional distribution, it is also possible to generate the starting values from the stationary distribution of a regime. See the details section.

a size [p, d, R2] array specifying the initial values in each slice for each Monte Carlo repetition, where d is the number of time series in the system and R2 is an argument of this function. In each slice, the last row will be used as initial values for the first lag, the second last row for second lag etc. If not specified, initial values will be drawn from the regime specified in init_regimes.

a numeric vector with values in 1,...,d (d=ncol(data)) specifying which the variables for which the impulse responses should be cumulative. Default is none.

should the GIRFs to some of the shocks be scaled so that they correspond to a specific magnitude of instantaneous or peak response of some specific variable (see the argument scale_type)? Provide a length three vector where the shock of interest is given in the first element (an integer in 1,...,d), the variable of interest is given in the second element (an integer in 1,...,d), and the magnitude of its instantaneous or peak response in the third element (a non-zero real number). If the GIRFs of multiple shocks should be scaled, provide a matrix which has one column for each of the shocks with the columns being the length three vectors described above.

If argument scale is specified, should the GIRFs be scaled to match an instantaneous response ("instant") or peak response ("peak"). If "peak", the scale is based on the largest magnitude of peak response in absolute value. Ignored if scale is not specified.

If scale_type == "peak" what the maximum horizon up to which peak response is expected? Scaling won't based on values after this.

a numeric vector with elements in (0, 1) specifying the confidence levels of the "confidence intervals" that do not quantify uncertainty about the true parameter value but only uncertainty about the initial value (and possibly sign and size of the shock) within the given regime.

set TRUE for a special feature in which for every possible length p history in the data, or a subset of them if so specified in the argument data_girf_pars, the GIRF is estimated for a shock that has the sign and size of the corresponding structural shock recovered from the data. If used, the argument which_shocks must specify only one shock. See the details section.

a length five numeric vector with the following elements determining settings for use_data_shocks=TRUE (concerns the single shock specified in the argument which_shocks):

1. An integer between 0 and M determining the (dominant) regime for which the GIRF should be calculated (0 for all regimes).

2. A number between 0.5 and 1 determining how large transition weight a regime should have to be considered dominant in a given time period (i.e., determining which histories are used to calculate the GIRF if the first element is not 0).

3. Either 0, -1, or 1, determining whether the GIRF should be calculated using shocks of all signs, only negative shocks, or only positive shocks, respectively.

4. Either, 0, 1, or 2, determining whether the GIRF should be calculated using shocks of all sizes, only small shocks, or only large shocks, respectively.

5. A strictly positive real number determining what size shocks are considered large and what size small "in the scale of standard deviations" (for example, if set to 2, shocks larger than that are considered large and shocks smaller than that are considered small; note that the standard deviations of the shocks are normalized to unity).

the number CPU cores to be used in parallel computing. Only single core computing is supported if an initial value is specified (and the GIRF won't thus be estimated multiple times).

Burn-in period for simulating initial values from a regime.

if weight_function="exogenous", provide a size (N+1 \times M) matrix of exogenous transition weights for the regimes: [h, m] for the (after-the-impact) period h-1 and regime m weight ([1, m] is for the impact period). Ignored if weight_function != "exogenous".

A numeric vector initializing the seeds for the random number generator for estimation of each GIRF. Should have the length of at least (extra seeds are removed from the end of the vector)...

If initial values are drawn using init_regime:

R2

If initial values are specified in init_values:

dim(init_values)[3]

If use_data_shocks=TRUE:

1 (the vector of seeds are generated according on the number of histories in the data that satisfy the conditions given in the argument data_girf_pars).

Set NULL for not initializing the seed.

employ parallel computing? If FALSE, does not print anything.

object of class 'girf' generated by the function GIRF.

numeric vector of length four that adjusts the [bottom_marginal, left_marginal, top_marginal, right_marginal] as the relative sizes of the marginals to the figures of the responses.

graphical parameters passed to plot method plotting the GIRFs

the number of decimals to print

an integer specifying the horizon how far to print the estimates and confidence intervals. The default is that all the values are printed.

an object of class 'stvar', created by, e.g., fitSTVAR or fitSSTVAR.

a numeric vector of length at most d (=ncol(data)) and elements in 1,...,d specifying the structural shocks for which the GIRF should be estimated.

a non-zero scalar value specifying the common size for all scalar components of the structural shock. Note that the conditional covariance matrix of the structural shock is normalized to an identity matrix and that the (generalized) impulse responses may not be symmetric with respect to the sign and size of the shock.

a positive integer specifying the horizon how far ahead should the generalized impulse responses be calculated.

the number of repetitions used to estimate GIRF for each initial value.

the number of initial values to use, i.e., to draw from init_regime if init_values are not specified. The confidence bounds will be sample quantiles of the GIRFs based on different initial values. Ignored if the argument init_value is specified. @param init_regime an integer in 1,...,M specifying the regime from which the initial values should be generated from (see ?simulate.stvar). If use_data_shocks=TRUE this is argument not used and data_girf_pars should be specified instead.

an integer in 1,...,M specifying the regime from which the initial values should be generated from (using a simulation procedure with a burn-in period). For models with Gaussian conditional distribution, it is also possible to generate the starting values from the stationary distribution of a regime. See the details section.

a size [p, d, R2] array specifying the initial values in each slice for each Monte Carlo repetition, where d is the number of time series in the system and R2 is an argument of this function. In each slice, the last row will be used as initial values for the first lag, the second last row for second lag etc. If not specified, initial values will be drawn from the regime specified in init_regimes.

a numeric vector with values in 1,...,d (d=ncol(data)) specifying which the variables for which the impulse responses should be cumulative. Default is none.

should the GIRFs to some of the shocks be scaled so that they correspond to a specific magnitude of instantaneous or peak response of some specific variable (see the argument scale_type)? Provide a length three vector where the shock of interest is given in the first element (an integer in 1,...,d), the variable of interest is given in the second element (an integer in 1,...,d), and the magnitude of its instantaneous or peak response in the third element (a non-zero real number). If the GIRFs of multiple shocks should be scaled, provide a matrix which has one column for each of the shocks with the columns being the length three vectors described above.

If argument scale is specified, should the GIRFs be scaled to match an instantaneous response ("instant") or peak response ("peak"). If "peak", the scale is based on the largest magnitude of peak response in absolute value. Ignored if scale is not specified.

If scale_type == "peak" what the maximum horizon up to which peak response is expected? Scaling won't based on values after this.

a numeric vector with elements in (0, 1) specifying the confidence levels of the "confidence intervals" that do not quantify uncertainty about the true parameter value but only uncertainty about the initial value (and possibly sign and size of the shock) within the given regime.

set TRUE for a special feature in which for every possible length p history in the data, or a subset of them if so specified in the argument data_girf_pars, the GIRF is estimated for a shock that has the sign and size of the corresponding structural shock recovered from the data. If used, the argument which_shocks must specify only one shock. See the details section.

a length five numeric vector with the following elements determining settings for use_data_shocks=TRUE (concerns the single shock specified in the argument which_shocks):

1. An integer between 0 and M determining the (dominant) regime for which the GIRF should be calculated (0 for all regimes).

2. A number between 0.5 and 1 determining how large transition weight a regime should have to be considered dominant in a given time period (i.e., determining which histories are used to calculate the GIRF if the first element is not 0).

3. Either 0, -1, or 1, determining whether the GIRF should be calculated using shocks of all signs, only negative shocks, or only positive shocks, respectively.

4. Either, 0, 1, or 2, determining whether the GIRF should be calculated using shocks of all sizes, only small shocks, or only large shocks, respectively.

5. A strictly positive real number determining what size shocks are considered large and what size small "in the scale of standard deviations" (for example, if set to 2, shocks larger than that are considered large and shocks smaller than that are considered small; note that the standard deviations of the shocks are normalized to unity).

the number CPU cores to be used in parallel computing. Only single core computing is supported if an initial value is specified (and the GIRF won't thus be estimated multiple times).

Burn-in period for simulating initial values from a regime.

if weight_function="exogenous", provide a size (N+1 \times M) matrix of exogenous transition weights for the regimes: [h, m] for the (after-the-impact) period h-1 and regime m weight ([1, m] is for the impact period). Ignored if weight_function != "exogenous".

A numeric vector initializing the seeds for the random number generator for estimation of each GIRF. Should have the length of at least (extra seeds are removed from the end of the vector)...

If initial values are drawn using init_regime:

R2

If initial values are specified in init_values:

dim(init_values)[3]

If use_data_shocks=TRUE:

1 (the vector of seeds are generated according on the number of histories in the data that satisfy the conditions given in the argument data_girf_pars).

Set NULL for not initializing the seed.

employ parallel computing? If FALSE, does not print anything.

a list parameters used for calculating counterfactual GIRFs (set to NULL for regular GIRFs). Contains the elements:

cfact_metatype

should be always "counterfactual_girf".

cfact_type

a character string indicating the type of counterfactual to be computed: should the path of the policy variable be fixed to some hypothetical path (cfact_type="fixed_path") in given impulse response horizons or should the responses of the policy variable to lagged and contemporaneous movements of some given variable be muted (cfact_type="muted_response")? See details for more information.

policy_var

a positive integer between 1 and d indicating the index of the policy variable considered in the counterfactual scenario. Note that policy_var is assumed to satisfy !(policy_var %in% which_shocks).

mute_var

a positive integer between 1 and d indicating the index of the variable to whose movements the policy variable specified in the argument policy_var should not react to in the counterfactual scenario. This indicates also the index of the shock to which the policy variable should not react to. It is assumed that mute_var != policy_var. This argument is only used when cfact_type="muted_response".

cfact_start

a positive integer between 1 and nsteps indicating the starting impulse response horizon period for the counterfactual behavior of the specified policy variable.

cfact_end

a positive integer between cfact_start and nsteps indicating the ending period for the counterfactual behavior of the specified policy variable.

cfact_path

a numeric vector of length cfact_end-cfact_start+1 indicating the hypothetical path of the policy variable specified in the argument policy_var. This argument is only used when cfact_type="fixed_path".

Important: If this is something else than NULL, it will change how the function behaves!

a (T \times d) matrix such that the ith row contains the vector y_{i}=(y_{1i},...,y_{di}) (dx1). That is, the initial values are excluded but the last observations is included.

a (T \times d) matrix such that the ith row contains the conditional mean of the process \mu_{y,i}.

a (d \times d \times M) array such that the slice [, , m] contains the conditional covariance matrix of regime m.

a (T \times M) matrix such that [t, m] contains the time t transition weights of the mth regime.

a (T \times dp) matrix such that the i:th row contains the vector (y_{i-1}',...,y_{i-p}') ((dp)x1), where y_{i}=(y_{1i},...,y_{di}) (dx1). That is, the initial values are included but the last observations not.

the ((dp)x1) mean vector, rep(all_mu[,m], times=p), that is the same for all observations.

the (dp \times dp) covariance matrix that is the same for all observations.

an object of class 'stvar' generated by fitSTVAR or STVAR, containing the freely estimated model.

an object of class 'stvar' generated by fitSTVAR or STVAR, containing the constrained model.

an object of class 'stvar' generated by fitSTVAR or STVAR.

a strictly positive integer specifying the number of lags to be tested.

should test for remaining autocorrelation or heteroskedasticity be calculated?

an object of class 'stvar' generated by fitSTVAR or STVAR, containing the model specified by the null hypothesis (i.e., the constrained model).

a matrix or class 'ts' object with d>1 columns. Each column is taken to represent a single times series. NA values are not supported. Ignore if defining a model without data is desired.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

number of times series in the system, i.e. ncol(data). This can be used to define STVAR models without data and can be ignored if data is provided.

a real valued vector specifying the parameter values. Should have the form \theta = (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu), where (see exceptions below):

• \phi_{m} = the (d \times 1) intercept (or mean) vector of the mth 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 mth 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},...,\phi_{M} 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} is the intercept parameter, A_{m,i} denotes the ith coefficient matrix of the mth regime, \Omega_{m} denotes the positive definite error term covariance matrix of the mth regime, and B_m is the invertible (d\times d) impact matrix of the mth regime. \nu_m is the degrees of freedom parameter of the mth regime. If parametrization=="mean", just replace each \phi_{m} 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.

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

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

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.

Perform penalized LS estimation that minimizes penalized RSS in which estimates close to breaking or not satisfying the usual stability condition are penalized? If TRUE, the tuning parameter is set by the argument penalty_params[2], and the penalization starts when the eigenvalues of the companion form AR matrix are larger than 1 - penalty_params[1].

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

If TRUE, estimates not satisfying the stability condition are allowed. Always FALSE if weight_function="relative_dens".

should approximate standard errors be calculated?

object of class 'stvar'.

currently not used.

number of digits to be printed.

if set to TRUE, instead of printing the estimates, prints the approximate standard errors using square roots of the diagonal of inverse of the observed information matrix.

an object of class 'stvar'.

should the series be plotted with the estimated transition weights or conditional means?

if set to TRUE then the print will include log-likelihood and information criteria values.

The input time series data.

A list describing the model structure.

The parameters of the model.

Approximate standard errors of the parameters, if calculated.

The transition weights of the model.

Conditional means of the regimes, if data is provided.

Total conditional means of the model, if data is provided.

Total conditional covariances of the model, if data is provided.

A list of unconditional moments including regime autocovariances, variances, and means.

Raw residuals, if data is provided.

Standardized residuals, if data is provided.

Recovered structural shocks, if applicable.

Log-likelihood of the model, if data is provided.

The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.

The parameter estimates from all estimation rounds, if applicable.

The log-likelihood of the estimates from all estimation rounds, if applicable.

Indicators of which estimation rounds converged, if applicable.

Indicators of which round of optimization each estimate belongs to, if applicable.

The seeds used in the estimation in fitSTVAR, if applicable.

The least squares estimates of the parameters in the form (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\alpha (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

a (T \times d) matrix such that the ith row contains the vector y_{i}=(y_{1i},...,y_{di}) (dx1). That is, the initial values are excluded but the last observations is included.

a (T \times d) matrix such that the ith row contains the conditional mean of the process \mu_{y,i}.

a (d \times d \times M) array such that the slice [, , m] contains the conditional covariance matrix of regime m.

a (T \times M) matrix such that [t, m] contains the time t transition weights of the mth regime.

the degrees of freedom parameter value (assumed larger than two).

a positive integer specifying the autoregressive order

the number of time series in the system, i.e., the dimension

[d, d, p] array containing the AR coefficient matrices

the (d\times d) positive definite error term covariance matrix

an object of class 'stvar' generated by fitSTVAR or STVAR, containing the model specified by the alternative hypothesis (i.e., the unconstrained model).

a size (k\timesn_params) matrix with full row rank specifying a part of the null hypothesis where n_params is the number of parameters in the (unconstrained) model. See details for more information.

a length k vector specifying a part of the null hypothesis. See details for more information.

a size (dxd) square matrix to be vectorized.

a vector containing the elements to be tested.

object of class "stvar"

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.

based on which estimation round should the model be constructed? An integer value in 1,...,nrounds. If specified, then which_largest is ignored.

should approximate standard errors be calculated?

The input time series data.

A list describing the model structure.

The parameters of the model.

Approximate standard errors of the parameters, if calculated.

The transition weights of the model.

Conditional means of the regimes, if data is provided.

Total conditional means of the model, if data is provided.

Total conditional covariances of the model, if data is provided.

A list of unconditional moments including regime autocovariances, variances, and means.

Raw residuals, if data is provided.

Standardized residuals, if data is provided.

Recovered structural shocks, if applicable.

Log-likelihood of the model, if data is provided.

The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.

The parameter estimates from all estimation rounds, if applicable.

The log-likelihood of the estimates from all estimation rounds, if applicable.

Indicators of which estimation rounds converged, if applicable.

Indicators of which round of optimization each estimate belongs to, if applicable.

The seeds used in the estimation in fitSTVAR, if applicable.

The least squares estimates of the parameters in the form (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\alpha (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

object of class "stvar"

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.

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.

the number of cores to be used in parallel computing.

logical: should the progress of the algorithm be printed?

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.

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.

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.

the number of cores to be used in parallel computing.

logical: should the progress of the algorithm be printed?

the degrees of freedom parameter value, a numeric scalar strictly larger than two.

the skewness parameter value, a numeric scalar strictly between -1 and 1.

a numeric vector specifying the point where the gradient or Hessian should be calculated.

a function that takes in argument x as the first argument.

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

object of class "stvar"

an object of class 'stvar' defining a structural or reduced form STVAR model. For a reduced form model (that is not readily identified statiscally), the shocks are automatically identified by the lower triangular Cholesky decomposition.

how many steps ahead should be predicted, i.e., the forecast horizon?

to how many independent simulations should the forecast be based on?

a numeric vector specifying the confidence levels of the prediction intervals.

should the pointforecast be based on sample "median" or "mean"?

if weight_function="exogenous", provide a size (nsteps x M) matrix of exogenous transition weights for the regimes: [step, m] for step steps ahead and regime m weight. Ignored if weight_function!="exogenous".

a character string indicating the type of counterfactual to be computed: should the path of the policy variable be fixed to some hypothetical path (cfact_type="fixed_path") in given forecast horizons or should the responses of the policy variable to lagged and contemporaneous movements of some given variable be muted (cfact_type="muted_response")? See details for more information.

a positive integer between 1 and d indicating the index of the policy variable considered in the counterfactual scenario.

a positive integer between 1 and d indicating the index of the variable to whose movements the policy variable specified in the argument policy_var should not react to in the counterfactual scenario. This indicates also the index of the shock to which the policy variable should not react to. It is assumed that mute_var != policy_var. This argument is only used when cfact_type="muted_response".

a positive integer between 1 and nsteps indicating the starting forecast horizon period for the counterfactual behavior of the specified policy variable.

a positive integer between cfact_start and nsteps indicating the ending period for the counterfactual behavior of the specified policy variable.

a numeric vector of length cfact_end-cfact_start+1 indicating the hypothetical path of the policy variable specified in the argument policy_var. This argument is only used when cfact_type="fixed_path".

object of class 'cfactfore' created by the function cfact_fore.

parameters passed to print.stvarpred printing the forecast.

a positive integer specifying the number of observations to be plotted along with the forecast.

should forecasts for transition weights be plotted?

how many significant digits to print?

an object of class 'stvar', created by, e.g., fitSTVAR or fitSSTVAR.

a numeric vector of length at most d (=ncol(data)) and elements in 1,...,d specifying the structural shocks for which the GIRF should be estimated.

a non-zero scalar value specifying the common size for all scalar components of the structural shock. Note that the conditional covariance matrix of the structural shock is normalized to an identity matrix and that the (generalized) impulse responses may not be symmetric with respect to the sign and size of the shock.

a positive integer specifying the horizon how far ahead should the generalized impulse responses be calculated.

the number of repetitions used to estimate GIRF for each initial value.

the number of initial values to use, i.e., to draw from init_regime if init_values are not specified. The confidence bounds will be sample quantiles of the GIRFs based on different initial values. Ignored if the argument init_value is specified. @param init_regime an integer in 1,...,M specifying the regime from which the initial values should be generated from (see ?simulate.stvar). If use_data_shocks=TRUE this is argument not used and data_girf_pars should be specified instead.

an integer in 1,...,M specifying the regime from which the initial values should be generated from (using a simulation procedure with a burn-in period). For models with Gaussian conditional distribution, it is also possible to generate the starting values from the stationary distribution of a regime. See the details section.

a size [p, d, R2] array specifying the initial values in each slice for each Monte Carlo repetition, where d is the number of time series in the system and R2 is an argument of this function. In each slice, the last row will be used as initial values for the first lag, the second last row for second lag etc. If not specified, initial values will be drawn from the regime specified in init_regimes.

a numeric vector with values in 1,...,d (d=ncol(data)) specifying which the variables for which the impulse responses should be cumulative. Default is none.

should the GIRFs to some of the shocks be scaled so that they correspond to a specific magnitude of instantaneous or peak response of some specific variable (see the argument scale_type)? Provide a length three vector where the shock of interest is given in the first element (an integer in 1,...,d), the variable of interest is given in the second element (an integer in 1,...,d), and the magnitude of its instantaneous or peak response in the third element (a non-zero real number). If the GIRFs of multiple shocks should be scaled, provide a matrix which has one column for each of the shocks with the columns being the length three vectors described above.

If argument scale is specified, should the GIRFs be scaled to match an instantaneous response ("instant") or peak response ("peak"). If "peak", the scale is based on the largest magnitude of peak response in absolute value. Ignored if scale is not specified.

If scale_type == "peak" what the maximum horizon up to which peak response is expected? Scaling won't based on values after this.

a numeric vector with elements in (0, 1) specifying the confidence levels of the "confidence intervals" that do not quantify uncertainty about the true parameter value but only uncertainty about the initial value (and possibly sign and size of the shock) within the given regime.

set TRUE for a special feature in which for every possible length p history in the data, or a subset of them if so specified in the argument data_girf_pars, the GIRF is estimated for a shock that has the sign and size of the corresponding structural shock recovered from the data. If used, the argument which_shocks must specify only one shock. See the details section.

a length five numeric vector with the following elements determining settings for use_data_shocks=TRUE (concerns the single shock specified in the argument which_shocks):

1. An integer between 0 and M determining the (dominant) regime for which the GIRF should be calculated (0 for all regimes).

2. A number between 0.5 and 1 determining how large transition weight a regime should have to be considered dominant in a given time period (i.e., determining which histories are used to calculate the GIRF if the first element is not 0).

3. Either 0, -1, or 1, determining whether the GIRF should be calculated using shocks of all signs, only negative shocks, or only positive shocks, respectively.

4. Either, 0, 1, or 2, determining whether the GIRF should be calculated using shocks of all sizes, only small shocks, or only large shocks, respectively.

5. A strictly positive real number determining what size shocks are considered large and what size small "in the scale of standard deviations" (for example, if set to 2, shocks larger than that are considered large and shocks smaller than that are considered small; note that the standard deviations of the shocks are normalized to unity).

the number CPU cores to be used in parallel computing. Only single core computing is supported if an initial value is specified (and the GIRF won't thus be estimated multiple times).

Burn-in period for simulating initial values from a regime.

if weight_function="exogenous", provide a size (N+1 \times M) matrix of exogenous transition weights for the regimes: [h, m] for the (after-the-impact) period h-1 and regime m weight ([1, m] is for the impact period). Ignored if weight_function != "exogenous".

A numeric vector initializing the seeds for the random number generator for estimation of each GIRF. Should have the length of at least (extra seeds are removed from the end of the vector)...

If initial values are drawn using init_regime:

R2

If initial values are specified in init_values:

dim(init_values)[3]

If use_data_shocks=TRUE:

1 (the vector of seeds are generated according on the number of histories in the data that satisfy the conditions given in the argument data_girf_pars).

Set NULL for not initializing the seed.

employ parallel computing? If FALSE, does not print anything.

a character string indicating the type of counterfactual to be computed: should the path of the policy variable be fixed to some hypothetical path (cfact_type="fixed_path") in given impulse response horizons or should the responses of the policy variable to lagged and contemporaneous movements of some given variable be muted (cfact_type="muted_response")? See details for more information.

a positive integer between 1 and d indicating the index of the policy variable considered in the counterfactual scenario. Note that policy_var is assumed to satisfy !(policy_var %in% which_shocks).

a positive integer between 1 and d indicating the index of the variable to whose movements the policy variable specified in the argument policy_var should not react to in the counterfactual scenario. This indicates also the index of the shock to which the policy variable should not react to. It is assumed that mute_var != policy_var. This argument is only used when cfact_type="muted_response".

a positive integer between 0 and N indicating the starting impulse response horizon period for the counterfactual behavior of the specified policy variable.

a positive integer between cfact_start and N indicating the ending period for the counterfactual behavior of the specified policy variable.

a numeric vector of length cfact_end-cfact_start+1 indicating the hypothetical path of the policy variable specified in the argument policy_var. This argument is only used when cfact_type="fixed_path".

object of class 'cfactgirf' created by the function cfact_girf.

parameters passed to print.stvargirf printing the girf.

how many significant digits to print?

an object of class 'stvar' defining a structural or reduced form STVAR model. For a reduced form model (that is not readily identified statiscally), the shocks are automatically identified by the lower triangular Cholesky decomposition.

a character string indicating the type of counterfactual to be computed: should the path of the policy variable be fixed to some hypothetical path (cfact_type="fixed_path") in given points of time or should the responses of the policy variable to lagged and contemporaneous movements of some given variable be muted (cfact_type="muted_response")? See details for more information.

a positive integer between 1 and d indicating the index of the policy variable considered in the counterfactual scenario.

a positive integer between 1 and d indicating the index of the variable to whose movements the policy variable specified in the argument policy_var should not react to in the counterfactual scenario. This indicates also the index of the shock to which the policy variable should not react to. It is assumed that mute_var != policy_var. This argument is only used when cfact_type="muted_response".

a positive integer between 1 and T indicating the starting period for the counterfactual behavior of the specified policy variable. It is assumed that the observed time series in indexed as y_{t-p+1},...,y_{0},y_1,...,y_T, i.e., that the first p observations are the initial values, and the "time period one" observation is the p+1th row in the data matrix.

a positive integer between cfact_start and T indicating the ending period for the counterfactual behavior of the specified policy variable.

a numeric vector of length cfact_end-cfact_start+1 indicating the hypothetical path of the policy variable specified in the argument policy_var. This argument is only used when cfact_type="fixed_path".

object of class 'cfacthist' created by the function cfact_hist.

arguments passed to the function window to select the time periods to print.

how many significant digits to print?

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

the number of time series in the system, i.e., the dimension

a real valued vector specifying the parameter values. Should have the form \theta = (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu), where (see exceptions below):

• \phi_{m} = the (d \times 1) intercept (or mean) vector of the mth 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 mth 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},...,\phi_{M} 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} is the intercept parameter, A_{m,i} denotes the ith coefficient matrix of the mth regime, \Omega_{m} denotes the positive definite error term covariance matrix of the mth regime, and B_m is the invertible (d\times d) impact matrix of the mth regime. \nu_m is the degrees of freedom parameter of the mth regime. If parametrization=="mean", just replace each \phi_{m} 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.

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

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

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

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

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

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.

If you want to switch between mean and intercept parametrizations:

either "intercept" or "mean" specifying to which parametrization it should be switched to. If set to "intercept", it's assumed that params is mean parametrized, and if set to "mean" it's assumed that params is intercept parametrized.

If you want to switch between the paramterizations B_{y,t}=\sum_{m=1}^M\alpha_{m,t}B_m and B_{y,t}=B_1 + \sum_{m=2}^M\alpha_{m,t}B_m^{*}, B_m^{*} = B_m - B_1:

either "orig" (with B_m) or "alt" (with B_m^{*}). It is assumed that the parameter vector is in the other parametrization than the one specified in change_to.

the autoregressive order of the model

the number of regimes

the number of time series in the system, i.e., the dimension

a real valued vector specifying the parameter values. Should have the form \theta = (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu), where (see exceptions below):

• \phi_{m} = the (d \times 1) intercept (or mean) vector of the mth 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 mth 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 threshold values.

weight_function="exogenous":

Omit \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.

Above, \phi_{m} is the intercept parameter, A_{m,i} denotes the ith coefficient matrix of the mth regime, \Omega_{m} denotes the positive definite error term covariance matrix of the mth regime, and B_m is the invertible (d\times d) impact matrix of the mth regime. \nu_m is the degrees of freedom parameter of the mth regime. If parametrization=="mean", just replace each \phi_{m} with regimewise mean \mu_{m}.

which regime?

If cond_dist="Gaussian" or cond_dist="Student":

rhe (dp + pd^2 + d(d+1)/2) vector (\phi_{m},vec(A_{m,1}),...,\vec(A_{m,p}),vech(\Omega_m)).

If cond_dist="ind_Student" or "ind_skewed_t":

the (dp + pd^2 + d^2 + 1) vector (\phi_{m},vec(A_{m,1}),...,\vec(A_{m,p}),vec(B_m)).

A cube (3D array) of impact matrices, with each slice being an inveritble square matrix.

A matrix of weights, with as many columns as there are slices in all_Omegas.

A strictly positive small number used as a tolerance for checking the invertibility of the matrix. The matrix is considered non-invertible if its determinant is less than this tolerance.

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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

the number of time series in the system, i.e., the dimension

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

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

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.

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.

a positive integer specifying the autoregressive order

Exogenous transition weights weights given for simulation in some context.

how many rows the exogenous weights should have?

what is the name of the object whose value should determine the the number of rows in the exogenous weights?

the autoregressive order of the model

the number of regimes

the number of time series in the system, i.e., the dimension

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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

the number of time series in the system, i.e., the dimension

a real valued vector specifying the parameter values. Should have the form \theta = (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu), where (see exceptions below):

• \phi_{m} = the (d \times 1) intercept (or mean) vector of the mth 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 mth 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},...,\phi_{M} 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} is the intercept parameter, A_{m,i} denotes the ith coefficient matrix of the mth regime, \Omega_{m} denotes the positive definite error term covariance matrix of the mth regime, and B_m is the invertible (d\times d) impact matrix of the mth regime. \nu_m is the degrees of freedom parameter of the mth regime. If parametrization=="mean", just replace each \phi_{m} 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.

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

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

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.

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

If TRUE, estimates not satisfying the stability condition are allowed. Always FALSE if weight_function="relative_dens".

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.

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.

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").

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.

S3 object to be tested

what is the name of the object that should of class 'stvar'?

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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

Integer, the lag for which the matrix is constructed, corresponding to the sub-diagonal filled with 1's.

Integer, the number of time periods or observations in the dataset, corresponding to the dimensions of the matrix.

An integer representing the dimension of the identity matrix.

An integer representing the factor by which to extend the matrix with zeros.

a positive definite (dxd) covariance matrix (d>1)

another positive definite (dxd) covariance matrix

object of class "stvar"

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.

should standardized or raw residuals be used?

the maximum lag considered in types "ac" and "ch".

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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

Perform penalized LS estimation that minimizes penalized RSS in which estimates close to breaking or not satisfying the usual stability condition are penalized? If TRUE, the tuning parameter is set by the argument penalty_params[2], and the penalization starts when the eigenvalues of the companion form AR matrix are larger than 1 - penalty_params[1].

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

the smallest accepted number of observations (times variables) from each regime relative to the number of parameters in the regime. For models with AR constraints, the number of AR matrix parameters in each regimes is simplisticly assumed to be ncol(AR_constraints)/M.

should the grid of weight function values in LS/NLS estimation be more sparse (speeding up the estimation)?

employ parallel computing? If FALSE, does not print anything.

the number CPU cores to be used in parallel computing.

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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

Perform penalized LS estimation that minimizes penalized RSS in which estimates close to breaking or not satisfying the usual stability condition are penalized? If TRUE, the tuning parameter is set by the argument penalty_params[2], and the penalization starts when the eigenvalues of the companion form AR matrix are larger than 1 - penalty_params[1].

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

the smallest accepted number of observations (times variables) from each regime relative to the number of parameters in the regime. For models with AR constraints, the number of AR matrix parameters in each regimes is simplisticly assumed to be ncol(AR_constraints)/M.

should the grid of weight function values in LS/NLS estimation be more sparse (speeding up the estimation)?

employ parallel computing? If FALSE, does not print anything.

the number CPU cores to be used in parallel computing.

a class 'stvar' object defining a structural STVAR model that is identified by heteroskedasticity or non-Gaussianity, typically created with fitSSTVAR.

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.

Should estimates close to breaking the usual stability condition be filtered out?

should approximate standard errors be calculated?

The input time series data.

A list describing the model structure.

The parameters of the model.

Approximate standard errors of the parameters, if calculated.

The transition weights of the model.

Conditional means of the regimes, if data is provided.

Total conditional means of the model, if data is provided.

Total conditional covariances of the model, if data is provided.

A list of unconditional moments including regime autocovariances, variances, and means.

Raw residuals, if data is provided.

Standardized residuals, if data is provided.

Recovered structural shocks, if applicable.

Log-likelihood of the model, if data is provided.

The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.

The parameter estimates from all estimation rounds, if applicable.

The log-likelihood of the estimates from all estimation rounds, if applicable.

Indicators of which estimation rounds converged, if applicable.

Indicators of which round of optimization each estimate belongs to, if applicable.

The seeds used in the estimation in fitSTVAR, if applicable.

The least squares estimates of the parameters in the form (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\alpha (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

a an object of class 'stvar', created by, e.g., fitSTVAR, specifying a reduced form or a structural model

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.

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.

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.

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

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

the maximum number of iterations in the variable metric algorithm.

the maximum number of iterations on the first phase robust estimation, if employed.

the strictly positive difference used in the finite difference approximation of the gradient used in numerical optimization.

Should some robust estimation method be used in the estimation before switching to the gradient based variable metric algorithm? See details.

should summaries of estimation results be printed?

Calculate approximate standard errors (based on standard asymptotics)?

The input time series data.

A list describing the model structure.

The parameters of the model.

Approximate standard errors of the parameters, if calculated.

The transition weights of the model.

Conditional means of the regimes, if data is provided.

Total conditional means of the model, if data is provided.

Total conditional covariances of the model, if data is provided.

A list of unconditional moments including regime autocovariances, variances, and means.

Raw residuals, if data is provided.

Standardized residuals, if data is provided.

Recovered structural shocks, if applicable.

Log-likelihood of the model, if data is provided.

The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.

The parameter estimates from all estimation rounds, if applicable.

The log-likelihood of the estimates from all estimation rounds, if applicable.

Indicators of which estimation rounds converged, if applicable.

Indicators of which round of optimization each estimate belongs to, if applicable.

The seeds used in the estimation in fitSTVAR, if applicable.

The least squares estimates of the parameters in the form (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\alpha (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

either "two-phase" or "three-phase" (the former is the default method for relative_dens models and the latter the default for the other models). See details.

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.

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

If TRUE, estimates not satisfying the stability condition are allowed. Always FALSE if weight_function="relative_dens".

In the LS/NLS step of the three phase estimation, the smallest accepted number of observations (times variables) from each regime relative to the number of parameters in the regime. For models with AR constraints, the number of AR matrix parameters in each regimes is simply assumed to be ncol(AR_constraints)/M.

should the grid of weight function values in LS/NLS estimation be more sparse (speeding up the estimation)?

the strictly positive difference used in the finite difference approximation of the gradient used in numerical optimization.

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

the number CPU cores to be used in parallel computing.

the maximum number of iterations in the variable metric algorithm.

a length nrounds vector containing the random number generator seed for each call to the genetic algorithm, or NULL for not initializing the seed.

should summaries of estimation results be printed?

employ parallel computing? If use_parallel=FALSE && print_res=FALSE, nothing is printed during the estimation process.

Calculate approximate standard errors (based on standard asymptotics)?

additional settings passed to the function GAfit employing the genetic algorithm.

The input time series data.

A list describing the model structure.

The parameters of the model.

Approximate standard errors of the parameters, if calculated.

The transition weights of the model.

Conditional means of the regimes, if data is provided.

Total conditional means of the model, if data is provided.

Total conditional covariances of the model, if data is provided.

A list of unconditional moments including regime autocovariances, variances, and means.

Raw residuals, if data is provided.

Standardized residuals, if data is provided.

Recovered structural shocks, if applicable.

Log-likelihood of the model, if data is provided.

The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.

The parameter estimates from all estimation rounds, if applicable.

The log-likelihood of the estimates from all estimation rounds, if applicable.

Indicators of which estimation rounds converged, if applicable.

Indicators of which round of optimization each estimate belongs to, if applicable.

The seeds used in the estimation in fitSTVAR, if applicable.

The least squares estimates of the parameters in the form (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\alpha (intercepts replaced by unconditional means if mean parametrization is used), if applicable.

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.

a positive integer specifying the autoregressive order

a positive integer specifying the number of regimes

a real valued vector specifying the parameter values. Should have the form \theta = (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\sigma,\alpha,\nu), where (see exceptions below):

• \phi_{m} = the (d \times 1) intercept (or mean) vector of the mth 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 mth 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},...,\phi_{M} 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} is the intercept parameter, A_{m,i} denotes the ith coefficient matrix of the mth regime, \Omega_{m} denotes the positive definite error term covariance matrix of the mth regime, and B_m is the invertible (d\times d) impact matrix of the mth regime. \nu_m is the degrees of freedom parameter of the mth regime. If parametrization=="mean", just replace each \phi_{m} 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.

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 mth 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 ith 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 ith 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 ith variable.

"exogenous":

Exogenous nonrandom transition weights, specify the weight series in weightfun_pars.

See the vignette for more details about the weight functions.

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.

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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters \phi_{m} or regime means \mu_{m}, m=1,...,M.

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

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

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

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

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.

A list containing internally used additional type of constraints (see the options below).

$fixed_lambdas (only if identification="heteroskedasticity"):

a length d(M-1) numeric vector (\lambda_{2},..., \lambda_{M}) with elements strictly larger than zero specifying the fixed parameter values for the parameters \lambda_{mi} should be constrained to.

$B1_constraints (only if identification="non-Gaussianity"):

set to the string "fixed_sign_and_order" to impose the constraints that the elements of the first impact matrix B_1 are strictly positive and that they are in a decreasing order.

Should some robust estimation method be used in the estimation before switching to the gradient based variable metric algorithm? See details.

Perform penalized LS estimation that minimizes penalized RSS in which estimates close to breaking or not satisfying the usual stability condition are penalized? If TRUE, the tuning parameter is set by the argument penalty_params[2], and the penalization starts when the eigenvalues of the companion form AR matrix are larger than 1 - penalty_params[1].

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

If TRUE, estimates not satisfying the stability condition are allowed. Always FALSE if weight_function="relative_dens".

the value that will be returned if the parameter vector does not lie in the parameter space (excluding the identification condition).

the maximum number of iterations in the variable metric algorithm.

the maximum number of iterations on the first phase robust estimation, if employed.

the seed for the random number generator (relevant when using SANN).

The input time series data.

A list describing the model structure.

The parameters of the model.

Approximate standard errors of the parameters, if calculated.

The transition weights of the model.

Conditional means of the regimes, if data is provided.

Total conditional means of the model, if data is provided.

Total conditional covariances of the model, if data is provided.

A list of unconditional moments including regime autocovariances, variances, and means.

Raw residuals, if data is provided.

Standardized residuals, if data is provided.

Recovered structural shocks, if applicable.

Log-likelihood of the model, if data is provided.

The values of the information criteria (AIC, HQIC, BIC) for the model, if data is provided.

The parameter estimates from all estimation rounds, if applicable.

The log-likelihood of the estimates from all estimation rounds, if applicable.

Indicators of which estimation rounds converged, if applicable.

Indicators of which round of optimization each estimate belongs to, if applicable.

The seeds used in the estimation in fitSTVAR, if applicable.

The least squares estimates of the parameters in the form (\phi_{1},...,\phi_{M},\varphi_1,...,\varphi_M,\alpha (intercepts replaced by unconditional means if mean parametrization is used), if applicable.