Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
55,421 characters · 7 sections · 0 citation commands
EM Estimation of Conditional Matrix Variate $t$ Distributions
An expectation--maximization (EM) algorithm was proposed and given its name by \citeA{Dempster77}. The EM algorithm is an iterative method to obtain (local) maximum likelihood estimates of parameters of distribution functions, which depend on unobserved (latent) variables. The EM algorithm alternates an expectation (E) step and a maximization (M) step. In the E--step, one considers that conditional on available data and the current estimate of the parameters, expectation of augmented log--likelihood of the data, and unobserved (latent) variables. The E--Step defines an objective function. In the M--step, to obtain a parameter estimate of the next iteration, one maximizes the objective function with respect to the parameters. Alternating between these steps, the EM algorithm produces improved parameter estimates at each step (in the sense that the value of the original log--likelihood is continually increased), and it converges to the maximum likelihood (ML) estimates of the parameters.
The EM algorithm is widely used in econometrics. In particular, \citeA{Hamilton90} introduced a parameter estimation method for a general regime-switching model. The regime-switching model assumes that a discrete unobservable Markov process randomly switches among a finite set of regimes and that a particular parameter set defines each regime. In finance, to value private companies whose market prices are unobservable, \citeA{Battulga23b} and \citeA{Battulga24c} applied the EM algorithm. \citeA{McNeil05} provides an EM algorithm to estimate parameters of the generalized hyperbolic distribution, which can be used to model financial returns. Also, the EM algorithm has been used in classifications. For example, to estimate matrix variate $t$ distribution parameters, \citeA{Thompson20} used the EM algorithm. \citeA{Sun10} provided an EM algorithm to estimate the parameters of Pearson VII distribution.
Classic Vector Autoregressive (VAR) process was proposed by \citeA{Sims80} who criticize large--scale macro--econometric models, which are designed to model interdependencies of economic variables. Besides \citeA{Sims80}, there are some other important works on multiple time series modeling, see, e.g., \citeA{Tiao81}, where a class of vector autoregressive moving average models was studied. For the VAR process, a variable in the process is modeled by its past values and the past values of other variables in the process. After the work of \citeA{Sims80}, VARs have been used for macroeconomic forecasting and policy analysis. However, if the number of variables in the system increases or the time lag is chosen high, then too many parameters need to be estimated. This will reduce the degrees of freedom of the model and entail a risk of over--parametrization.
Therefore, to reduce the number of parameters in a high--dimensional VAR process, \citeA{Litterman79} introduced probability distributions for coefficients that are centered at the desired restrictions but that have a small and nonzero variance. Those probability distributions are known as Minnesota prior in Bayesian VAR (BVAR) literature, which is widely used in practice. Due to over--parametrization, the generally accepted result is that the forecast of the BVAR model is better than the VAR model estimated by the frequentist technique. The BVAR relies on Monte--Carlo simulation. Recently, for Bayesian Markov--Switching VAR process, \citeA{Battulga24g} introduced a new Monte--Carlo simulation method that removes duplication in a regime vector. Also, the author introduced importance sampling method to estimate probability of rare event, which corresponds to endogenous variables. Research works have shown that BVAR is an appropriate tool for modeling large data sets; for example, see \citeA{Banbura10}.
The rest of the paper is organized as follows: In Section 2, for type I conditional matrix variate distributions, including general case and special case with Minnesota prior we develop EM algorithms. Section 3 is dedicated to studying EM algorithms for type II conditional matrix variate distributions, including general case and special case with Minnesota prior. Finally, Section 4 concludes the study.
We consider a Bayesian Vector Autoregressive process of $p$ order (BVAR($p$)), which is given by the following equation
where $y_t=(y_{1,t},\dots,y_{n,t})'$ is an $(n\times 1)$ vector of endogenous variables, $\psi_t=(1,\psi_{2,t},\dots,\psi_{l,t})'$ is an $(l\times 1)$ vector of exogenous variables, $\xi_t=(\xi_{1,t},\dots,\xi_{n,t})'$ is an $(n\times 1)$ residual process, $A_{0,t}$ is an $(n\times l)$ random coefficient matrix, corresponding to the vector of exogenous variables, for $i=1,\dots,p$, $A_{i,t}$ are $(n\times n)$ random coefficient matrices, corresponding to $y_{t-1},\dots,y_{t-p}$. Equation (ref) can be written by
where $\Pi_t:=[A_{0,t}: A_{1,t}:\dots:A_{p,t}]$ is an $(n\times d)$ random coefficient matrix with $d:=l+np$, which consist of all the random coefficient matrices and $\mathsf{Y}_t:=(\psi_t',y_{t-1}',\dots,y_{t-p}')'$ is a $(d\times 1)$ vector, which consist of exogenous variable $\psi_t$ and last $p$ lagged values of the process $y_t$. The process $\mathsf{Y}_t$ is measurable with respect to a $\sigma$--field $\mathcal{F}_{t-1}$, which is defined below. Let us collect endogenous variables into ($[nT]\times 1$) vector $y$, i.e., $y:=(y_1',\dots,y_T')'$.
For the residual process $\xi_t$, we assume that it has $\xi_t:=\Sigma_t^{1/2}\varepsilon_t$, $t=1,\dots,T$ representation, where $\Sigma_t^{1/2}$ is a Cholesky factor of a positive definite $(n\times n)$ random matrix $\Sigma_t$ and $\varepsilon_1,\dots,\varepsilon_T$ is a random sequence of independent identically multivariate normally distributed random vectors with means of 0 and covariance matrices of $n$ dimensional identity matrix $I_n$. We also assume that the strong white noise process $\{\varepsilon_t\}_{t=1}^T$ is independent of the random coefficient matrices $(\Pi_1,\dots,\Pi_T)$ and $(\Sigma_1,\dots,\Sigma_T)$ conditional on initial information $\mathcal{F}_0:=\{y_{1-p},\dots,y_0,\psi_{1},\dots,\psi_T\}$, where $\psi_1,\dots,\psi_T$ are values of exogenous variables and they are known at time zero. We also denote available information at time $t$ by $\mathcal{F}_t:=\{\mathcal{F}_0,y_1,\dots,y_t\}$.
Let us assume that for $t=1,\dots,T$, the random coefficient matrices $\Pi_t$ and random covariance matrices $\Sigma_t$ are equals, that is, $\Pi:=\Pi_1=\dots=\Pi_T$ and $\Sigma:=\Sigma_1=\dots=\Sigma_T$. Then, by using the Kronecker product, the BVAR($p$) process can be written by the following equation
where $\otimes$ is the Kronecker product of two matrices and $\pi:=\text{vec}(\Pi)$ is an $(nd\times 1)$ vectorization of the random coefficient matrix $\Pi$. Now we define distributions of the random coefficient vector $\pi$ and covariance matrix $\Sigma$. We assume that conditional on the initial information $\mathcal{F}_0$, a distribution of the random covariance matrix $\Sigma$ is given by
where the notation $\mathcal{IW}$ denotes the Inverse--Wishart distribution, $\nu_0>n-1$ is a degrees of freedom and $V_0$ is a positive definite scale matrix. Consequently, a distribution of the residual vector $\xi_t$ equals
where $\mathcal{N}$ denotes the normal distribution. Also, we assume that conditional on the covariance matrix $\Sigma$ and initial information $\mathcal{F}_0$, a distribution of the random coefficient vector $\pi$ is given by
where $\pi_0$ is an $(nd\times 1)$ vector and $\Lambda_{0}=\text{diag}\{\lambda_{1},\dots,\lambda_{d}\}$ is a diagonal $(d\times d)$ matrix. Then, according to \citeA{Battulga24g}, the following Proposition holds.
The joint density function (ref) is called conditional matrix variate student $t$ density, see \citeA{Battulga24g}. To differentiate from a conditional matrix $t$ distribution, which arise in the following section, we refer to the distribution as type I conditional matrix variate student $t$ distribution. From equations (ref) and (ref), one recognizes that posterior distributions of the random coefficient vector $\pi$ and covariance matrix $\Sigma$ are multivariate normal $\mathcal{N}\big(\pi_{0|t},\Lambda_{0|T}\otimes \Sigma\big)$ and inverse--Wishart $\mathcal{IW}\big(\nu_0+T,B_T+V_0\big)$, respectively. Let us denote a vector of all the parameters of the joint density function by $\theta:=\text{vec}(\pi_0,\Lambda_0,\nu_0,V_0)$. In this section, we develop Expectation--Maximization (EM) algorithm to estimate parameters of the density function. In E--Step, we consider that conditional on the full information $\mathcal{F}_T$ and parameter at iteration $k$, $\theta^{[k]}$, expectation of augmented log--likelihood of the data $y$ and unobserved (latent) variables $\pi$ and $\Sigma$. The E--Step defines a objective function $\mathcal{L}$, namely,
In M--Step, to obtain parameter estimate of next iteration $\theta^{[k+1]}$, one maximizes the objective function with respect to the parameter $\theta$. First, let us consider partial derivative from the objective function with respect to the parameter $\lambda_i$ for $i=1,\dots,d$. Since an inverse of the matrix $\Lambda_0$ is $\Lambda_0^{-1}=\text{diag}\{1/\lambda_1,\dots,1/\lambda_d\}$, it is clear that
where $\mathsf{E}_{ij}^d$ is a $(d\times d)$ matrix and its $(i,j)$--th element equals 1 and others 0 and for generic square matrix $A$, $\text{tr}(A)$ denotes trace of the matrix $A$. In general case, to obtain parameter estimation of the matrix $\Lambda_0$, one may use the following partial derivative
It follows from mean and covariance matrix of the random vector $\pi$, given in equation (ref) that
where
is a positive semi--definite $(d\times d)$ matrix,
is a positive semi--definite $(d\times d)$ matrix and
is an $([nd]\times 1)$ vector. $\Lambda_{0|T}^{[k]}$ and $\pi_{0|T}^{[k]}$ are Bayesian estimators at iteration $k$ of the parameter matrix $\Lambda_0$ and the parameter vector $\pi$ at iteration $k$, respectively. Note that the vector $\pi_{0|T}^{[k]}$ does not depend on the random covariance matrix $\Sigma$. According to the iterated expectation formula and expectation formula of Wishart distributed random matrix $\Sigma^{-1}$, we have that
where the matrix $B_T^{[k]}$ equals
It can be shown that the matrix $\Psi_{i|T}^{[k]}$ is a positive semi--definite matrix. Consequently, as $\ln|\Lambda_0|=\sum_{i=1}^d\ln(\lambda_i)$, for $i=1,\dots,d$, an estimator at iteration $(k+1)$ of the parameter $\lambda_i$ is given by
Because $\Psi_{i|T}^{[k]}$ is a positive semi--definite matrix, $\lambda_i^{[k+1]}$ takes a non--negative value. For $i=1,\dots,d$, we collect $\lambda_i^{[k+1]}$ into a matrix $\Lambda_0^{[k+1]}:=\text{diag}\Big\{\lambda_1^{[k+1]},\dots,\lambda_d^{[k+1]}\Big\}$.
Second, we consider partial derivative from the objective function $\mathcal{L}$ with respect to the parameter vector $\pi_0$. The partial derivative is given by
According to the iterated expectation formula and equations (ref) and (ref), one obtains estimator at iteration $(k+1)$ of the parameter vector $\pi_0$
Third, we consider partial derivative from the objective function $\mathcal{L}$ with respect to the parameter matrix $V_0$. The partial derivative is given by
Consequently, we have that
Fourth, we consider partial derivative from the objective function $\mathcal{L}$ with respect to the parameter $\nu_0$. The partial derivative is given by
where $\psi_n(\cdot)$ is the multivariate digamma function, which is defined by derivative of log of the multivariate gamma function. As a result, since $\mathbb{E}\big[\ln|\Sigma^{-1}|\big|\mathcal{F}_T\big]=\psi_n(\frac{\nu_0+T}{2})+n\ln(2)-\ln|V_0+B_T|$, see \citeA{Nguyen23}, one gets that
If we substitute equation (ref) into the above equation, then we have that
As a result, to obtain an estimator at iteration $(k+1)$ of the parameter $\nu_0$, one has to solve the above nonlinear equation for $\nu_0$.
In the following algorithm, we give EM algorithm for parameters of conditional matrix variate $t$ density function
In this subsection, we follow \citeA{Battulga24g} and we develop EM algorithm for parameters of type I conditional matrix variate $t$ distribution with Minnesota prior. In practice, one usually adopts the Minnesota prior to estimating the parameters of the VAR$(p)$ process. The first version of Minnesota prior was introduced by \citeA{Litterman79}. Also, \citeA{Banbura10} used Minnesota prior for large Bayesian VAR and showed that the forecast of large Bayesian VAR is better than small Bayesian VAR. However, there are many different variants of the Minnesota prior, we consider a prior, which is included in \citeA{Miranda18}. The idea of Minnesota prior is that it shrinks diagonal elements of the matrix $A_1$ toward $\phi_i$ and off--diagonal elements of $A_1$ and all elements of other matrices $A_0,A_2,\dots,A_p$ toward 0, where $\phi_i$ is 0 for a stationary variable $y_{i,t}$ and 1 for a variable with unit root $y_{i,t}$. However, we adopt a different prior condition for the random coefficient matrix $A_0$. Without loss of generality let us assume that there are $m$ $(m=0,\dots,n)$ stationary variables and the stationary variables are placed on the first $m$ components of the process $y_t$. For the prior, it is assumed that conditional on $\Sigma$ and $\mathcal{F}_0$, $A_0,A_1,\dots,A_p$ are jointly normally distributed, and for $(i,j)$--th element of the matrix $A_{\ell}$ $(\ell=0,\dots,p)$, it holds that for $i=1,\dots,n$ and $j=1,\dots,l$,
and
and for $i,j=1,\dots,n$,
and
where $C_{i,j}$ is $(i,j)$--th element of an $(n\times m)$ parameter matrix $C_m$ and $\sigma_i^2$ is an $(i,i)$--th element of the random covariance matrix $\Sigma$. It follows from equation (ref) that the expectation of $(A_0)_{i,j}$ equals $C_{i,j}$ for stationary variable $y_{i,t}$ and zero for non stationary variable with unit root $y_{i,t}$. Let us introduce an $(n\times l)$ matrix $\bar{C}_m:=[C_m:0_{[n\times (l-m)]}]$. A small $\varepsilon_i^2$ corresponds to an uninformative diffuse prior for $(A_0)_{i,j}$, the parameter $\alpha$ controls the overall tightness of the prior distribution, the parameter $\beta$ controls amount of information prior information at higher lags, and $\gamma_i$ is a scaling parameter, see \citeA{Miranda18}. Thus, the factor $1/\ell^{2\beta}$ represents a rate at which prior variance decreases with increasing lag length.
According to \citeA{Banbura10}, it can be shown that the following equation satisfies the prior conditions (ref)--(ref)
where $\hat{y}^\circ$ and $\hat{\mathsf{Y}}^\circ$ are $(n\times d)$ and $(d\times d)$ matrices of dummy variables and are defined by
and
with $J_\beta:=\text{diag}\{1^{\beta},\dots,p^{\beta}\}$, respectively, and $\hat{\xi}^\circ:=[\xi_1:\dots:\xi_d]$ is an $(n\times d)$ matrix of residual process. Note that one can add constraints for elements of the coefficient matrix $\Pi$ to the matrices of dummy variables. It is worth mentioning that the matrices of dummy variables $\hat{y}_t$ and $\hat{\mathsf{Y}}_t$ should not depend on the covariance matrix $\Sigma$. If the dummy variables depend on the covariance matrix, an OLS estimator, and matrix $\Lambda_{0}$ depend on the covariance matrix $\Sigma$, see below. Consequently, in this case, one can not use the results of Proposition (ref). For this reason, we choose the prior condition (ref)--(ref). Equation (ref) can be written by
where $\hat{y}$ and $\hat{\xi}$ are $([nd]\times 1)$ vectors and are vectorizations of the matrix of dummy variables $\hat{y}^\circ$ and matrix of the residual process $\hat{\xi}^\circ$, respectively, i.e., $\hat{y}:=\text{vec}(\hat{y}^\circ)$ and $\hat{\xi}:=\text{vec}(\hat{\xi}^\circ)$. It follows from equation (ref) that
where $d$ denotes equal distribution. It should be noted that the first term of the right--hand side of the above equation is a vecorization of the ordinary least square (OLS) estimator of the coefficient matrix $\Pi$, namely,
Note that the OLS estimator of the coefficient matrix $\Pi$ are same as the prior conditions (ref) and (ref). Consequently, conditional on $\Sigma$ and $\mathcal{F}_0$, a distribution of the coefficient vector $\pi$ is given by
where $\Lambda_0:=(\hat{\mathsf{Y}}^\circ(\hat{\mathsf{Y}}^\circ)')^{-1}$ is a $(d\times d)$ diagonal matrix and its inverse equals
Since the matrix $\Lambda_0^{-1}$ is a diagonal matrix, its determinant, which appears in the objective function (ref) is
Let for $m=1,\dots,l$, $c_m:=\text{vec}(C_m)$ be a $(nm\times 1)$ vectorization of the parameter matrix $C_m$. Then, partial derivative from the objective function $\mathcal{L}$ with respect to the parameter $c$ is
where $J_m:=\big[I_{mn}:0_{[mn\times n^2p]}\big]$ is an $(mn\times dn)$ matrix and $D_m(\varepsilon):=\text{diag}\{\varepsilon_1^2,\dots,\varepsilon_m^2\}$ is an $(m\times m)$ diagonal matrix. Then, since $J_m(\pi-\pi_0)=J_m\pi-c_m$, we get an estimator at iteration $(k+1)$ of the parameter vector $c_m$
where $\pi_{0|T}^{[k]}$ is given by equation (ref). Note that if all variables of the process $y_t$ are unit root processes ($m=0$), then one does not need a parameter estimation of the matrix $C_m$. To obtain estimators at iteration $(k+1)$ of the parameters $\alpha$, $\beta$, $\gamma_i$, and $\varepsilon_j$ for $i=1,\dots,n$ and $j=1,\dots,l$, we define the following $(d\times d)$ matrices:
and
It follows from determinant equation (ref) and the objective function (ref) that similarly to equation (ref), one obtains that for $i=1,\dots,n$ and $j=1,\dots,l$, estimators at iteration $(k+1)$ of the parameters $\varepsilon_j$, $\alpha$, and $\gamma_i$ are given by
and
where $\Theta_{|T}^{[k]}$ is calculated via equations (ref), (ref), and (ref). For $p\geq 2$, an estimator at iteration $(k+1)$ of the parameter $\beta$, $\beta^{[k+1]}$ is obtained from the following nonlinear equation
It should be noted that it is not difficult to show that the right--hand sides of equations (ref)--(ref) take non negative values. Consequently, an EM algorithm for parameters of conditional matrix variate $t$ distribution with Minnesota prior is given by the following algorithm.
Let us reconsider a BVAR($p$) process, given in equation (ref), namely,
where $\pi_t:=\text{vec}(\Pi_t)$ is an $(nd\times 1)$ vectorization of the random coefficient matrix $\Pi_t$. For the random coefficient vectors and covariance matrices, we assume the following assumption holds.
Note that from the assumption, we can conclude that conditional on initial information $\mathcal{F}_0$, for each $t=2,\dots,T$, $(\pi_t,\Sigma_t)$ is independent of a random vector $(y_1',\dots,y_{t-1}')'$. Now we define distributions of the random coefficient vector $\pi_t$ and covariance matrix $\Sigma_t$. We suppose that conditional on the initial information $\mathcal{F}_0$, a distribution of the random covariance matrix $\Sigma_t$ is given by
where $\nu_0>n-1$ is a degrees of freedom and $V_0$ is a positive definite scale matrix. Hence, a distribution of the residual vector $\xi_t$ equals
Also, we suppose that conditional on the covariance matrix $\Sigma_t$ and initial information $\mathcal{F}_0$, a distribution of the random coefficient vector $\pi_t$ is given by
where $\pi_0$ is an $(nd\times 1)$ vector and $\Lambda_0$ is a symmetric positive definite $(d\times d)$ matrix. Then, the following Proposition holds.
We refer to distribution function, which is given in equation (ref) as type II conditional matrix variate $t$ distribution function. In this section, we introduce an EM algorithm to estimate the parameters of the distribution function. Similarly to equation (ref), an conditional expectation (objective function) in E--Step is given by
In M--Step, to obtain ML estimators at iteration $(k+1)$ of the parameters of the type II conditional matrix variate $t$ distribution, we need the following Lemma.
By using Proposition (ref), Lemma (ref), and ideas in subsection 2.1, one can arrive the following EM Algorithm, which estimates parameters of the general type II conditional matrix variate $t$ distribution.
Now we consider type II conditional matrix variate $t$ distribution with Minnesota prior. For random coefficient matrices $A_{0,t},A_{1,t},\dots,A_{p,t}$, we assume that same prior conditions hold as subsection 2.2, namely, for $i=1,\dots,n$ and $j=1,\dots,l$,
and
and for $i,j=1,\dots,n$,
and
where $\sigma_{i,t}^2$ is an $(i,i)$--th element of the random covariance matrix $\Sigma_t$. Then, similarly to the Algorithm (ref), one obtains the following EM algorithm, which estimate parameters of type II conditional matrix variate $t$ distribution with Minnesota prior.
Conditional matrix variate student $t$ distribution was introduced by \citeA{Battulga24g}. In this paper, we provide EM algorithms, which estimate parameters of the conditional matrix variate student $t$ distributions, including general case and special case with Minnesota prior. Also, we introduce a new conditional matrix variate student $t$ distribution, which is closely related to the \citeA{Battulga24g}'s conditional matrix variate student $t$ distribution.
Here we provide proofs of Proposition (ref) and Lemma (ref).