fitGMVAR: Two-phase maximum likelihood estimation of a GMVAR model

Description

fitGMVAR estimates a GMVAR 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

fitGMVAR(
  data,
  p,
  M,
  conditional = TRUE,
  parametrization = c("intercept", "mean"),
  constraints = NULL,
  structural_pars = NULL,
  ncalls = round(10 + 9 * log(M)),
  ncores = min(2, ncalls, parallel::detectCores()),
  maxit = 300,
  seeds = NULL,
  print_res = TRUE,
  ...
)

Arguments

data

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

a positive integer specifying the autoregressive order of the model.

a positive integer specifying the number of mixture components.

conditional

a logical argument specifying whether the conditional or exact log-likelihood function

parametrization

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

constraints

a size $(Mpd^2 x q)$ constraint matrix $C$ specifying general linear constraints to the autoregressive parameters. We consider constraints of form ($\phi$$_{1}$$,...,$$\phi$$_{M}) = $$C \psi$, where $\phi$$_{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 same for all regimes, set $C$= [I:...:I]' $(Mpd^2 x pd^2)$ where I = diag(p*d^2). Ignore (or set to NULL) if linear constraints should not be employed.

structural_pars

If NULL a reduced form model is considered. For structural model, should be a list containing the following elements:

W - a $(dxd)$ matrix with its entries imposing constraints on $W$: 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.
C_lambda - a $(d(M-1) x r)$ constraint matrix that satisfies ($\lambda$$_{2}$$,...,$ $\lambda$$_{M}) =$ $C_{\lambda} \gamma$ where $\gamma$ is the new $(r x 1)$ parameter subject to which the model is estimated (similarly to AR parameter constraints). The entries of C_lambda must be either positive or zero. Ignore (or set to NULL) if the eigenvalues $\lambda_{mi}$ should not be constrained.

See Virolainen (2020) for the conditions required to identify the shocks and for the B-matrix as well (it is $W$ times a time-varying diagonal matrix with positive diagonal entries).

ncalls

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 ncalls vector containing the random number generator seed for each call to the genetic algorithm, or NULL for not initializing the seed. Exists for creating reproducible results.

print_res

should summaries of estimation results be printed?

...

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

Value

Returns an object of class 'gmvar' defining the estimated (reduced form or structural) GMVAR model. Multivariate quantile residuals (Kalliovirta and Saikkonen 2010) are also computed and included in the returned object. In addition, the returned object contains the estimates and log-likelihood values from all the estimation rounds performed. The estimated parameter vector can be obtained at gmvar$params (and corresponding approximate standard errors at gmvar$std_errors) and it is...

For unconstrained models:

...a size $((M(pd^2+d+d(d+1)/2+1)-1)x1)$ vector that has form $\theta$$ = $($\upsilon$$_{1}$, ...,$\upsilon$$_{M}$, $\alpha_{1},...,\alpha_{M-1}$), where

$\upsilon$$_{m}$ $ = (\phi_{m,0},$$\phi$$_{m}$$,\sigma_{m})$
$\phi$$_{m}$$ = (vec(A_{m,1}),...,vec(A_{m,p})$
and $\sigma_{m} = vech(\Omega_{m})$, m=1,...,M.

For constrained models:

...a size $((M(d+d(d+1)/2+1)+q-1)x1)$ vector that has form $\theta$$ = (\phi_{1,0},...,\phi_{M,0},$$\psi$ $,\sigma_{1},...,\sigma_{M},\alpha_{1},...,\alpha_{M-1})$, where

$\psi$ $(qx1)$ satisfies ($\phi$$_{1}$$,...,$ $\phi$$_{M}) =$ $C \psi$ where $C$ is $(Mpd^2xq)$ constraint matrix.

For structural GMVAR model:

...a vector that has the form $\theta$$ = (\phi_{1,0},...,\phi_{M,0},$$\phi$$_{1},...,$$\phi$$_{M}, vec(W),$$\lambda$$_{2},...,$$\lambda$$_{M},\alpha_{1},...,\alpha_{M-1})$, where

$\lambda$$_{m}=(\lambda_{m1},...,\lambda_{md})$ contains the eigenvalues of the $m$th mixture component.

If AR parameters are constrained:: Replace $\phi$$_{1}$$,...,$ $\phi$$_{M}$ with $\psi$ $(qx1)$ that satisfies ($\phi$$_{1}$$,...,$ $\phi$$_{M}) =$ $C \psi$, as above.
If $W$ is constrained:: Remove the zeros from $vec(W)$ and make sure the other entries satisfy the sign constraints.
If $\lambda_{mi}$ are constrained:: Replace $\lambda$$_{2},...,$$\lambda$$_{M}$ with $\gamma$ $(rx1)$ that satisfies ($\lambda$$_{2}$$,...,$ $\lambda$$_{M}) =$ $C_{\lambda} \gamma$ where $C_{\lambda}$ is a $(d(M-1) x r)$ constraint matrix.

Above, $\phi_{m,0}$ is the intercept parameter, $A_{m,i}$ denotes the $i$th coefficient matrix of the $m$th mixture component, $\Omega_{m}$ denotes the error term covariance matrix of the $m$:th mixture component, and $\alpha_{m}$ is the mixing weight parameter. The $W$ and $\lambda_{mi}$ are structural parameters replacing the error term covariance matrices (see Virolainen, 2020). If $M=1$, $\alpha_{m}$ and $\lambda_{mi}$ are dropped. If parametrization=="mean", just replace each $\phi_{m,0}$ with regimewise mean $\mu_{m}$. $vec()$ is vectorization operator that stacks columns of a given matrix into a vector. $vech()$ stacks columns of a given matrix from the principal diagonal downwards (including elements on the diagonal) into a vector. The notation is in line with the cited article by Kalliovirta, Meitz and Saikkonen (2016) introducing the GMVAR model.

Remark that the first autocovariance/correlation matrix in $uncond_moments is for the lag zero, the second one for the lag one, etc.

S3 methods

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

Details

Because of complexity and multimodality of the log-likelihood function, it's not certain that the estimation algorithms will end up in the global maximum point. It's expected that most of the estimation rounds will end up in some local maximum or saddle point instead. Therefore, a (sometimes large) number of estimation rounds is required for reliable results. Because of the nature of the model, the estimation may fail especially in the cases where the number of mixture components is chosen too large.

The estimation process is computationally heavy and it might take considerably long time for large models with large number of observations. 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 and employ another round of estimation (see ?GAfit how to set up an initial population with the dot parameters).

The code of the genetic algorithm is mostly based on the description by Dorsey and Mayer (1995) but it includes some extra features that were found useful for this particular estimation problem. For instance, the genetic algorithm uses a slightly modified version of the individually adaptive crossover and mutation rates described by Patnaik and Srinivas (1994) and employs (50%) fitness inheritance discussed by Smith, Dike and Stegmann (1995).

The gradient based variable metric algorithm used in the second phase is implemented with function optim from the package stats.

Finally, note that the structural models are even more difficult to estimate than the reduced form models due to the different parametrization of the covariance matrices. If necessary, an initial population may be constructed for the genetic algorithm based on the estimation results of a reduced form model. It does not seem unambiguous, however, how to do that so that the (structural) parameter constraints are satisfied. Also, be aware that if the lambda parameters are constrained in any other way than by restricting some of them to be identical, the parameter "lambda_scale" of the genetic algorithm (see ?GAfit) needs to be carefully adjusted accordingly. Be aware that if a structural model is considered and the lambda parameters are constrained in some other way than constraining some of them to be identical, the settings of the genetic algorithm may have to be adjusted.

References

Dorsey R. E. and Mayer W. J. 1995. Genetic algorithms for estimation problems with multiple optima, nondifferentiability, and other irregular features. Journal of Business & Economic Statistics, 13, 53-66.
Kalliovirta L., Meitz M. and Saikkonen P. 2016. Gaussian mixture vector autoregression. Journal of Econometrics, 192, 485-498.
Kalliovirta L. and Saikkonen P. 2010. Reliable Residuals for Multivariate Nonlinear Time Series Models. Unpublished Revision of HECER Discussion Paper No. 247.
Patnaik L.M. and Srinivas M. 1994. Adaptive Probabilities of Crossover and Mutation in Genetic Algorithms. Transactions on Systems, Man and Cybernetics 24, 656-667.
Smith R.E., Dike B.A., Stegmann S.A. 1995. Fitness inheritance in genetic algorithms. Proceedings of the 1995 ACM Symposium on Applied Computing, 345-350.
Virolainen S. 2020. Structural Gaussian mixture vector autoregressive model. Unpublished working paper, available as arXiv:2007.04713.

Examples

Run this code

# NOT RUN {
## These are long running examples that use parallel computing!

# These examples use the data 'eurusd' which comes with the
# package, but in a scaled form (similar to Kalliovirta et al. 2016).
data(eurusd, package="gmvarkit")
data <- cbind(10*eurusd[,1], 100*eurusd[,2])
colnames(data) <- colnames(eurusd)

# GMVAR(1,2) model: 10 estimation rounds with seeds set
# for reproducibility
fit12 <- fitGMVAR(data, p=1, M=2, ncalls=10, seeds=1:10)
fit12
plot(fit12)
summary(fit12)
print_std_errors(fit12)
profile_logliks(fit12)

# Structural GMVAR(1,2) model identified with sign
# constraints. The sign constraints (which fully identify
# the shocks) are in line with the reduced form model,
# so the maximized loglikelihood is the same.
W_122 <- matrix(c(1, NA, -1, 1), nrow=2)
fit12s <- fitGMVAR(data, p=1, M=2, structural_pars=list(W=W_122),
  ncalls=10, seeds=1:10)
fit12s

# GMVAR(2,2) model with mean parametrization
fit22 <- fitGMVAR(data, p=2, M=2, parametrization="mean",
                  ncalls=16, seeds=1:16)
fit22

# Structural GMVAR(2,2) model with the lambda parameters restricted
# to be identical (in the second regime) and the shocks identified
# with diagonal of the B-matrix normalized positive and one zero constraint.
# The resulting model has error term covariance matrices that are
# multiplicatives of each other, while the identification equals to
# identification through Cholesky decomposition.
W_222 <- matrix(c(1, NA, 0, 1), nrow=2)
C_lambda_222 <- matrix(c(1, 1), nrow=2)
fit22s <- fitGMVAR(data, p=2, M=2, structural_pars=list(W=W_222, C_lambda=C_lambda_222),
  ncalls=20, seeds=1:20)
fit22s

# GMVAR(2,2) model with autoregressive parameters restricted
# to be the same for both regimes
C_mat <- rbind(diag(2*2^2), diag(2*2^2))
fit22c <- fitGMVAR(data, p=2, M=2, constraints=C_mat)
fit22c

# GMVAR(2,2) model with autoregressive parameters restricted
# to be the same for both regimes and non-diagonl elements
# the coefficient matrices constrained to zero. Estimation
# with only 10 estimation rounds.
tmp <- matrix(c(1, rep(0, 10), 1, rep(0, 8), 1, rep(0, 10), 1),
 nrow=2*2^2, byrow=FALSE)
C_mat2 <- rbind(tmp, tmp)
fit22c2 <- fitGMVAR(data, p=2, M=2, constraints=C_mat2)
fit22c2
# }

Run the code above in your browser using DataLab