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,
  same_means = NULL,
  structural_pars = NULL,
  ncalls = floor(10 + 30 * log(M)),
  ncores = min(2, ncalls, parallel::detectCores()),
  maxit = 500,
  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

"intercept" or "mean" determining whether the model is parametrized with intercept parameters $\phi_{m,0}$ or regime means $\mu_{m}$, 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.

same_means

Restrict the mean parameters of some regimes to be the same? 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 the same 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".

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

If you wish to estimate a structural model without overidentifying constraints that is identified statistically, specify your W matrix is structural_pars to be such that it contains the same sign constraints in a single row (e.g. a row of ones) and leave the other elements as NA. In this way, the genetic algorithm works the best. The ordering and signs of the columns of the W matrix can be changed afterwards with the functions reorder_W_columns and swap_W_signs.

Because of complexity and high 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).

If the estimation algorithm fails to create an initial population for the genetic algorithm, it usually helps to scale the individual series so that the AR coefficients (of a VAR model) will be relative small, preferably less than one. Even if one is able to create an initial population, it should be preferred to scale the series so that most of the AR coefficients will not be very large, as the estimation algorithm works better with 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. If initial population is still not found, you can try to adjust the parameters of the genetic algorithm according to the characteristics of the time series (for the list of the available settings, see ?GAfit).

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.

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, so larger number of estimation rounds should be considered. 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.

Finally, the function fails to calculate approximative standard errors and the parameter estimates are near the border of the parameter space, it might help to use smaller numerical tolerance for the stationarity and positive definiteness conditions. The numerical tolerance of an existing model can be changed with the function update_numtols.

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.
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!
# Running all the below examples will take approximately 3-4 minutes.

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

# The rest of the examples only use a single estimation round with a given
# seed that produces the MLE to reduce running time of the examples. When
# estimating models for empirical applications, a large number of estimation
# rounds (ncalls = a large number) should be performed to ensure reliability
# of the estimates.

# Structural GMVAR(1,2) model identified with sign
# constraints.
W_122 <- matrix(c(1, 1, -1, 1), nrow=2)
fit12s <- fitGMVAR(data, p=1, M=2, structural_pars=list(W=W_122),
  ncalls=1, seeds=1)
fit12s

# Structural GMVAR(2, 2) model identified statistically only
W_222 <- matrix(c(1, NA, 1, NA), nrow=2)
fit22s <- fitGMVAR(data, p=2, M=2, structural_pars=list(W=W_222),
  ncalls=1, seeds=12)
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, ncalls=1, seeds=1)
fit22c

# GMVAR(2,2) model with autoregressive parameters restricted
# to be the same for both regimes and non-diagonal 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,
                    ncalls=1, seeds=3)
fit22c2
# }

Run the code above in your browser using DataLab