Introduction to gmvarkit

Savi Virolainen

2019-08-28

The package gmvarkit contains tools to estimate and analyze Gaussian Mixture Vector Autoregressive (GMVAR) model (see Kalliovirta et. al. 2016). We briefly give the definition of the GMVAR model in this vignette, but we do not cover its (very attractive) theoretical and practical properties. The reader is hence adviced to familiarize oneself with the cited article(s).

The GMVAR models in gmvarkit are defined as class gmvar S3 objects whose can be created with the estimation function fitGMVAR or with the constructor function GMVAR. The created gmvar objects are then be conveniently used as main arguments in many other functions that allow, for example, model diagnostics, simulations and forecasting. Therefore, after estimating a GMVAR model it’s easy to use the other functions for further analyses. However some tasks, such as applying general linear constraints, creating GMVAR models without estimation or setting up initial population for the genetic algorithm employed by the estimation function, require accurate understanding on how the parameter vectors are constructed in this package.

The notations for (unconstrained) parameter vector are in line with the cited article by Kalliovirta et. al. (2016), but for clarity we repeat these notations in the second section of this vignette. In the third section we show how to apply general linear constraints to the autoregressive parameters of a GMVAR model. Two examples are given to demonstrate how to use the general formula. In the fourth section we have listed some useful functions and methods found in gmvarkit and briefly explain how they work and how to use them.

Definition of the GMVAR model

The Gaussian mixture vector autoregressive (GMVAR) model with \(M\) mixture components and autoregressive order \(p\), which we refer to as the GMVAR(\(p,M\)) model, is a mixture of \(M\) stationary Gaussian VAR(\(p\)) models with time varying mixing weights. At each point of time \(t\), a GMVAR process generates an observation from one of its mixture components (which are linear VAR models) according to probabilities pointed by the mixing weights.

Let \(y_t=(y_{1t},...,y_{dt})\), \(t=1,2,...\) be the \(d\)-dimensional vector valued time series of interest, and let \(\mathcal{F}_{t-1}\) be the \(\sigma\)-algebra generated by the random variables \(\lbrace y_{t-j},j>0 \rbrace\) (containing the information on the past of \(y_t\)). Let \(M\) be the number of mixture components, and \(\boldsymbol{s}_t=(s_{t,1},...,s_{t,M})\) a sequence of (non-observable) random vectors such that at each \(t\) exactly one of its components takes the value one and others take the value zero. The component that takes the value one is selected according to the probabilities \(Pr(s_{t,m}=1|\mathcal{F}_{t-1})\equiv\alpha_{m,t}\), \(m=1,...,M\) with \(\sum_{m=1}^M\alpha_{m,t}=1\) for all \(t\). Then, for a GMVAR(\(p,M\)) model, we have \[ y_t = \sum_{m=1}^M s_{t,m}(\mu_{m,t}+\Omega_{m}^{1/2}\varepsilon_t), \quad \varepsilon_t\sim\text{NID}(0,I_d), \] where \(\Omega_{m}\) is a conditional positive definite covariance matrix, NID stands for “normally and identically distributed”, and \[ \mu_{m,t} = \varphi_{m,0} + \sum_{i=1}^p A_{m,i}y_{t-i} \] is interpreted as the conditional mean of the \(m\)th mixture components. The terms \(\varphi_{m,0}\) \((d\times 1)\) are intercept parameters and \(A_{m,i}\) \((d\times d)\) are the autoregressive (AR) matrices. The random vectors \(\varepsilon_t\) are assumed to be independent from \(\mathcal{F}_{t-1}\) and conditionally indepedent from \(\boldsymbol{s}_{t}\) given \(\mathcal{F}_{t-1}\).

For each \(m=1,...,M\), the AR matrices \(A_{m,i}\), \(i=1,...,p\) are assumed to satisfy the usual stability condition of VAR models. Namely, denoting \[ \boldsymbol{A}_m= \begin{bmatrix} A_1 & A_2 & \cdots & A_{p-1} & A_p \\ I_d & 0 & \cdots & 0 & 0 \\ 0 & I_d & \cdots & 0 & 0 \\ \vdots & & \ddots & 0 & 0 \\ 0 & 0 & \cdots & I_d & 0 \end{bmatrix} \enspace (dp\times dp) \] we assume that all the modulus of the eigenvalues of the matrix \(\boldsymbol{A}_m\) are smaller than one, for each \(m=1,...,M\).

Based on the above definition, the conditional density function of \(y_t\) given \(\mathcal{F}_{t-1}\), \(f(y_t|\mathcal{F}_{t-1})\), is given as \[ f(y_t|\mathcal{F}_{t-1}) = \sum_{m=1}^M\alpha_{m,t}(2\pi)^{-d/2}\det(\Omega_m)^{-1/2}\exp\left\lbrace -\frac{1}{2}(y_t-\mu_{m,t})´\Omega_m^{-1}(y_t-\mu_{m,t}) \right\rbrace . \] In order to define the mixing weights \(\alpha_{m,t}\), consider first auxiliary stationary Gaussian VAR-processes \(z_{m,t}\), and the related \(dp\)-dimensional random vectors \(\boldsymbol{z}_{t,m}=(z_{t,m},...,z_{m,t-p+1})\). The density of \(\boldsymbol{z}_{t,m}\) is \[ n_{dp}(\boldsymbol{z}_{t,m};\mu_m,\boldsymbol{\Sigma}_{m,p})=(2\pi)^{-dp/2}\det(\boldsymbol{\Sigma}_{m,p})^{-1/2}\exp\left\lbrace -\frac{1}{2}(\boldsymbol{z}_{t,m} - \boldsymbol{1}_p\otimes\mu_m)´\boldsymbol{\Sigma}_{m,p}^{-1}(\boldsymbol{z}_{t,m} - \boldsymbol{1}_p\otimes\mu_m)\right\rbrace , \] where \(\boldsymbol{1}_p = (1,...,1)\) \((p\times 1)\), \(\mu_m=(I_d - \sum_{i=1}^pA_{m,i})^{-1}\varphi_{m,0}\), and \(\boldsymbol{\Sigma}_{m,p}\) is obtained from \(vec(\boldsymbol{\Sigma}_{m,p})=(I_{dp^2} - \boldsymbol{A}_m\otimes\boldsymbol{A}_m)^{-1}vec(\Sigma_{m,\varepsilon})\) where \[ \Sigma_{m,\varepsilon} = \begin{bmatrix} \Omega_m & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \\ \end{bmatrix} \enspace (dp\times dp). \]

Denoting \(\boldsymbol{y}_{t-1}=(y_{t-1},...,y_{t-p})\) \((dp\times 1)\), the mixing weights of the GMVAR model are defined as \[ \alpha_{m,t} = \frac{\alpha_mn_{dp}(\boldsymbol{y}_{t-1};\mu_m,\boldsymbol{\Sigma}_{m,p})}{\sum_{n=1}^M\alpha_nn_{dp}(\boldsymbol{y}_{t-1};\mu_n,\boldsymbol{\Sigma}_{n,p})}, \] where the mixing weights parameters \(\alpha_m\) are assumed to satisfy \(\sum_{m=1}^M\alpha_{m}=1\). The probabilities for each regime occuring therefore depend on the level, variability, and temporal depedence of the past \(p\) observations.

General notations and unconstrained models

General notations

Defining a GMVAR model requires specifying the autoregressive degree p and the number of mixture components M. The number of time series in the system is denoted by d, and it’s assumed that d is larger than one.1 The form of the parameter vector depends on whether an unconstrained or constrained model is considered.

Parameter vector of unconstrained model

The parameter vector for unconstrained model is of size ((M(pd^2+d(d+1)/2+1)-1)x1) and has form \[\boldsymbol{\theta} = (\boldsymbol{\upsilon_1},...,\boldsymbol{\upsilon_M},\alpha_1,...,\alpha_{M-1}),\enspace \text{where}\quad\quad\enspace\quad\] \[\boldsymbol{\upsilon_m}=(\phi_{m,0},\boldsymbol{\phi_m},\sigma_m)\enspace ((pd^2+d(d+1)/2)x1),\quad\quad\enspace\] \[\boldsymbol{\phi_m}=(vec(A_{m,1}),...,vec(A_{m,p}))\enspace (pd^2x1), \enspace \text{and}\quad\enspace\;\] \[\sigma_m=vech(\Omega_m)\enspace ((d(d+1)/2)x1),\enspace m=1,...,M.\]

Above \(\phi_{m,0}\) denotes the intercept parameter of m:th mixture component, \(A_{m,i}\) denotes the coefficient matrix of m:th mixture component and i:th lag, \(\Omega_m\) denotes the (positive definite) error term covariance matrix, and \(\alpha_m\) denotes the mixing weight parameter of m:th mixture component. \(vec()\) is a vectorization operator that stacks columns of a matrix into a vector and \(vech()\) stacks columns of a matrix from the main diagonal downwards (including the main diagonal) into a vector.

The parameter vector above has “intercept” parametrization, referring to the intercept terms \(\phi_{m,0}\). However, it’s also possible to use “mean” parametrization, where the intercept terms are simply replaced by the regimewise means \(\mu_m=(I_d-\sum_{i=1}^pA_{m,i})^{-1}\phi_{m,0}\).

Models with linear constraints

General linear constraints and parameter vector

Imposing linear constraints on the autoregressive parameters of GMVAR model is straightforward in gmvarkit. The constraints are expressed in a rather general form allowing arbitrary linear constraints, but one needs to take the time to construct the constraint matrix carefully for each particular case.

We consider constraints of form \[(\boldsymbol{\phi_1},...,\boldsymbol{\phi_M}) = \boldsymbol{C}\boldsymbol{\psi},\enspace \text{where}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\enspace\] \[\boldsymbol{\phi_m}=(vec(A_{m,1}),...,vec(A_{m,p}))\enspace (pd^2x1), \enspace m=1,...,M,\] \(\boldsymbol{C}\) is known \((Mpd^2xq)\) constraint matrix (of full column rank) and \(\boldsymbol{\psi}\) is unknown \((qx1)\) parameter vector.

The parameter vector for constrained model is size ((M(d+d(d+1)/2+1)+q-1)x1) and has form \[\boldsymbol{\theta} = (\phi_{1,0},...,\phi_{M,0},\boldsymbol{\psi},\alpha_1,...,\alpha_{M-1}),\] where \(\boldsymbol{\psi}\) is the \((qx1)\) parameter vector containing constrained autoregressive parameters. As in the case of regular models, instead of the intercept parametrization that takes use of intercept terms \(\phi_{m,0}\), one may use the mean parametrization with regimewise means \(\mu_m\) instead (m=1,…,M).

Examples of linear constraints

Consider the following two common uses for linear constraints: restricting the autoregressive parameters to be the same for all regimes and constraining some parameters to zero. Of course also some other constraints may be useful, but we chose to show illustrative examples of these two as they occur in the cited article by Kalliovirta et al. (2016).

Restricting AR parameters to be the same for all regimes

To restrict the AR parameters to be the same for all regimes, we want \(\boldsymbol{\phi_m}\) to be the same for all m=1,…,M. The parameter vector \(\boldsymbol{\psi}\) \((qx1)\) then corresponds to any \(\boldsymbol{\phi_m}=\boldsymbol{\phi}\), and therefore \(q=pd^2\). For the constraint matrix we choose \[\boldsymbol{C} = [I_{pd^2}:\cdots:I_{pd^2}]' \enspace (Mpd^2xpd^2),\] that is, M pieces of \((pd^2xpd^2)\) diagonal matrices stacked on top of each other, because then \[\boldsymbol{C}\boldsymbol{\psi}=(\boldsymbol{\psi},...,\boldsymbol{\psi})=(\boldsymbol{\phi},...,\boldsymbol{\phi}).\]

Restricting AR parameters to be the same for all regimes and constraining non-diagonal elements of coefficient matrices to be zero

The previous example shows how to restrict the AR parameters to be the same for all regimes, but say we also want to constrain the non-diagonal elements of coefficient matrices \(A_{m,i}\) (m=1,…,M, i=1,…,p) to be zero. We have the constrained parameter \(\boldsymbol{\psi}\) \((qx1)\) representing the unconstrained parameters \((\boldsymbol{\phi_1},...,\boldsymbol{\phi_M})\), where by assumption \(\boldsymbol{\phi_m}=\boldsymbol{\phi}=(vec(A_1),...,vec(A_p))\) \((pd^2x1)\) and the elements of \(vec(A_i)\) (i=1,…,p) corresponding to the diagonal are zero.

For illustrative purposes, let’s consider a GMVAR model with autoregressive degree p=2, number of mixture components M=2 and number of time series in the system d=2. Then we have \[\boldsymbol{\phi}=(A_{1,(1,1)},0,0,A_{1,(2,2)},A_{2,(1,1)},0,0,A_{2,(2,2)}) \enspace (8x1) \enspace \text{and} \] \[\boldsymbol{\psi}=(A_{1,(1,1)},A_{1,(2,2)},A_{2,(1,1)},A_{2,(2,2)}) \enspace (4x1).\quad\quad\quad\quad\quad\enspace\] By a direct calculation, we can see that choosing the constraint matrix \[\boldsymbol{C}=\left[{\begin{array}{c} \boldsymbol{\tilde{c}} \\ \boldsymbol{\tilde{c}} \\ \end{array}}\right] \enspace (Mpd^2x4), \enspace \text{where}\]

\[\boldsymbol{\tilde{c}}=\left[{\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array}}\right] \enspace (pd^2x4)\] satisfies \(\boldsymbol{C}\boldsymbol{\psi}=(\boldsymbol{\phi},...,\boldsymbol{\phi}).\)

Functions in gmvarkit

Estimating a GMVAR model

gmvarkit provides the function fitGMVAR for estimating a GMVAR model. The maximum likelihood estimation is performed in two phases: in the first phase fitGMVAR uses a genetic algorithm to find starting values for gradient based variable metric algorithm, which it then uses to finalize the estimation in the second phase. It’s important to keep in mind that it’s not guaranteed that the numerical estimation algorithms end up in the global maximum point rather than a local one. Because of multimodality and challenging surface of the log-likelihood function, it’s actually expected that most of the estimation rounds won’t find the global maximum point. For this reason one should always perform multiple estimation rounds. Parallel computing is used to obtain the results faster. The number of CPU cores used can be set with the argument ncores, and the number of estimation rounds can be controlled with the argument ncalls.

If the model estimates poorly, it is often because the number of mixture components M is chosen too large. One may also adjust the settings of the genetic algorithm employed, or set up an initial population with guesses for the estimates. This can by done by passing arguments in fitGMVAR to the function GAfit which implements the genetic algorithm. To check the available settings, read the documentation ?GAfit. If the iteration limit is reached when estimating the model, the function iterate_more can be used to finish the estimation.

Parameters of the estimated model are printed in an illustrative and easy to read form. In order to easily compare approximate standard errors to certain estimates, one can print the approximate standard errors of the estimates in the same form with the function print_std_errors. Numerical approximation of the gradient and Hessian matrix of the log-likelihood at the estimates can be obtained conveniently with the functions get_gradient and get_hessian, and eigenvalues of the Hessian can be obtained with the function get_soc.

The estimated objects have their own print, plot, and summary methods.

The function alt_gmvar can be used to build GMVAR model based on an arbitrary estimation round. This can be used to consider estimates pointed by local maximums instead of the ML estimate.

Model diagnostics

gmvarkit considers model diagnostics based on multivariate extension of quantile residuals (see Kalliovirta and Saikkonen 2010), whose are under the correct model specification asymptotically multivariate standard normal distributed. The quantile residual tests introduced by Kalliovirta and Saikkonen (2010) can be performed with the function quantile_residual_tests by providing the estimated model (that is class gmvar object) as an argument.

For graphical diagnostics one may use the function diagostic_plot, which enables one to plot the quantile residual time series, auto- and cross-correlation functions of quantile residuals or their squares, or quantile residual histograms and normal QQ-plots.

Constructing a GMVAR model without estimation

One may wish to construct an arbitrary GMVAR model without any estimation process, for example in order to simulate from a particular process of interest. An arbitrary model can be created with the function GMVAR. If one wants to add or update data to the model afterwards, it’s advisable to use the function add_data.

Simulating from GMVAR process

The function simulateGMVAR is the one for the job. As the main argument it uses a gmvar object created with fitGMVAR or GMVAR.

Forecasting GMVAR process

The package gmvarkit contains predict method predict.gmvar for forecasting GMVAR processes. For one step predictions using the exact formula for conditional mean is supported, but the forecasts further than that are based on independent simulations. The predictions are either sample means or medians and the confidence intervals are based on sample quantiles. The objects generated by predict.gmvar have their own plot method. We also encourage directly using the function simulateGMVAR for forecasting.

Univariate analysis

Use the package uGMAR for analysing univariate time series.

References


  1. For univariate analysis one may use the package uGMAR.