fitSTVAR {sstvars}R Documentation

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

Description

fitSTVAR estimates a reduced form smooth transition VAR model in two phases: in the first phase, it uses a genetic algorithm to find starting values for a gradient based variable metric algorithm, which it then uses to finalize the estimation in the second phase. Parallel computing is utilized to perform multiple rounds of estimations in parallel.

Usage

fitSTVAR(
  data,
  p,
  M,
  weight_function = c("relative_dens", "logistic", "mlogit", "exponential", "threshold",
    "exogenous"),
  weightfun_pars = NULL,
  cond_dist = c("Gaussian", "Student", "ind_Student"),
  parametrization = c("intercept", "mean"),
  AR_constraints = NULL,
  mean_constraints = NULL,
  weight_constraints = NULL,
  nrounds = (M + 1)^5,
  ncores = 2,
  maxit = 1000,
  seeds = NULL,
  print_res = TRUE,
  use_parallel = TRUE,
  filter_estimates = TRUE,
  ...
)

Arguments

data

a matrix or class 'ts' object with d>1 columns. Each column is taken to represent a univariate time series. Missing values are not supported.

p

a positive integer specifying the autoregressive order

M

a positive integer specifying the number of regimes

weight_function

What type of transition weights \alpha_{m,t} should be used?

"relative_dens":

\alpha_{m,t}= \frac{\alpha_mf_{m,dp}(y_{t-1},...,y_{t-p+1})}{\sum_{n=1}^M\alpha_nf_{n,dp}(y_{t-1},...,y_{t-p+1})}, where \alpha_m\in (0,1) are weight parameters that satisfy \sum_{m=1}^M\alpha_m=1 and f_{m,dp}(\cdot) is the dp-dimensional stationary density of the 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.

weightfun_pars
If weight_function == "relative_dens":

Not used.

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

a numeric vector with the switching variable i\in\lbrace 1,...,d \rbrace in the first and the lag j\in\lbrace 1,...,p \rbrace in the second element.

If weight_function == "mlogit":

a list of two elements:

The first element $vars:

a numeric vector containing the variables that should used as switching variables in the weight function in an increasing order, i.e., a vector with unique elements in \lbrace 1,...,d \rbrace.

The second element $lags:

an integer in \lbrace 1,...,p \rbrace specifying the number of lags to be used in the weight function.

If weight_function == "exogenous":

a size (nrow(data) - p x M) matrix containing the exogenous transition weights as [t, m] for time t and regime m. Each row needs to sum to one and only weakly positive values are allowed.

cond_dist

specifies the conditional distribution of the model as "Gaussian", "Student", or "ind_Student", where the latest is the Student's t distribution with independent components.

parametrization

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

AR_constraints

a size (Mpd^2 x 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 x 1),\ m=1,...,M, contains the coefficient matrices and \psi (q x 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 x pd^2) where I = diag(p*d^2).

mean_constraints

Restrict the mean parameters of some regimes to be identical? Provide a list of numeric vectors such that each numeric vector contains the regimes that should share the common mean parameters. For instance, if M=3, the argument list(1, 2:3) restricts the mean parameters of the second and third regime to be identical but the first regime has freely estimated (unconditional) mean. Ignore or set to NULL if mean parameters should not be restricted to be the same among any regimes. This constraint is available only for mean parametrized models; that is, when parametrization="mean".

weight_constraints

a list of two elements, R in the first element and r in the second element, specifying linear constraints on the transition weight parameters \alpha. The constraints are of the form \alpha = R\xi + r, where R is a known (a\times l) constraint matrix of full column rank (a is the dimension of \alpha), r is a known (a\times 1) constant, and \xi is an unknown (l\times 1) parameter. Alternatively, set R=0 in order to constrain the the weight parameter to the constant r (in this case, \alpha is dropped from the constrained parameter vector).

nrounds

the number of estimation rounds that should be performed.

ncores

the number CPU cores to be used in parallel computing.

maxit

the maximum number of iterations in the variable metric algorithm.

seeds

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

print_res

should summaries of estimation results be printed?

use_parallel

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

filter_estimates

should the likely inappropriate estimates be filtered? See details.

...

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

Details

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

Because of complexity and high multimodality of the log-likelihood function, it is not certain that the estimation algorithm will end up in the global maximum point. It is expected that many of the estimation rounds will end up in some local maximum or a saddle point instead. Therefore, a (sometimes very large) number of estimation rounds is required for reliable results. Due to identification problems and high complexity of the surface of the log-likelihood function, the estimation may fail especially in the cases where the number of regimes is chosen too large.

The estimation process is computationally heavy and it might take considerably long time for large models to estimate. Note that estimation of model with cond_dist == "ind_Student" is computationally substantially more demanding than standard Gaussian and Student's t models due to the different structure of the log-likelihood function (parametrized directly via impact matrices rather than covariance matrices of the regimes).

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

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

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

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

Value

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

data

The input time series data.

model

A list describing the model structure.

params

The parameters of the model.

std_errors

Approximate standard errors of the parameters, if calculated.

transition_weights

The transition weights of the model.

regime_cmeans

Conditional means of the regimes, if data is provided.

total_cmeans

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

total_ccovs

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

uncond_moments

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

residuals_raw

Raw residuals, if data is provided.

residuals_std

Standardized residuals, if data is provided.

structural_shocks

Recovered structural shocks, if applicable.

loglik

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

IC

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

all_estimates

The parameter estimates from all estimation rounds, if applicable.

all_logliks

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

which_converged

Indicators of which estimation rounds converged, if applicable.

which_round

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

S3 methods

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

References

See Also

fitSSTVAR, STVAR, GAfit, iterate_more

Examples


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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