Variational Bayes Estimation of Discrete-Margined Copula Models with Application to Time Series
Ruben Loaiza-Maya Michael Stanley Smith
stat.ME econ.EM stat.ML


empty






Variational Bayes Estimation of Discrete-Margined Copula Models with Application to  Time Series
    Rubén Loaiza-Maya and Michael Stanley Smith
    First Version December 2017
This Version December 30, 2023
================================================================================================

Rubén Loaiza-Maya is a PhD student 
and Michael Smith is Chair of Management (Econometrics), both
at Melbourne Business School, University of Melbourne.


Correspondence should be directed to Michael Smith at
mike.smith@mbs.edu. 




Variational Bayes Estimation of Discrete-Margined Copula Models with Application to  Time Series



Abstract




We propose a new variational Bayes estimator for high-dimensional
copulas with discrete, or a combination of discrete and continuous, margins. The method
is based on a variational approximation to a tractable augmented posterior, and is faster than previous likelihood-based approaches.






We use it to estimate drawable vine copulas for univariate and
multivariate Markov ordinal and mixed time series.
These have dimension rT, where T is the number of observations
and r is the number of series,
and are difficult to estimate using previous methods. 
The vine pair-copulas are carefully selected to allow for heteroskedasticity, 
which is a feature of most
ordinal time series data. When combined with flexible margins, the resulting 
time series models also allow for other common features of ordinal data, such as
zero inflation, multiple modes and under- or over-dispersion.
Using six example series, we illustrate both the flexibility of the time series copula models, and  
the efficacy of the variational Bayes estimator for copulas of up to 792 dimensions 
and 60 parameters. 
This far exceeds the size and complexity of copula models for 
discrete data that can be estimated
using previous methods. 





Key Words: 
Data Augmentation; Drawable vines;  Heteroskedasticity; Multivariate ordinal and mixed time series;
Sparse variational approximation; Stochastic gradient ascent.


	




plain







§ INTRODUCTION



Copula models for m discrete-valued variables are difficult to estimate because their 
likelihood involves 2^m evaluations of the copula function, so that it is computationally intractable for even moderate dimensions. To avoid this problem, 
<cit.> and <cit.> propose using Markov chain Monte Carlo (MCMC)  data augmentation, where a tractable 
augmented likelihood is employed instead. 
However, this approach becomes slow
for copulas with higher dimensions and sample sizes.
<cit.> were the first to
suggest using a variational Bayes (VB) estimator as an alternative.
Their method — based on that in <cit.> and labeled VBIL — uses an unbiased estimate to the intractable likelihood computed
using importance sampling. However, as the copula dimension
or the number of copula parameters increase, computing the unbiased estimate of the likelihood
via importance sampling
also makes this method computationally infeasible. 

In this paper, we propose a new
VB estimator for copulas with a substantially
higher dimension and number of parameters
than can be estimated by either MCMC data augmentation or VBIL. It
uses a variational approximation to the tractable augmented likelihood of <cit.>,
instead of the intractable likelihood. We label our method
VBDA for `variational Bayes data augmentation', and propose several
variational approximations which  balance computational efficiency against accuracy. 
While the new VBDA method is a general
approach to estimate large copula models with one or more discrete margins, 
we employ it here to estimate 
multivariate times series copula models for ordinal, or a combination
of continuous and ordinal (ie. `mixed'), time series variables. The models combine
arbitrary
time-invariant margins with a copula that
captures serial and cross-sectional dependence jointly <cit.>. These
copulas are  
challenging to estimate because they have dimension rT, where T is the number of 
time series observations
and r is the number of series. 
We first show that VBDA is accurate when compared to the
exact posterior computed using
(much slower) MCMC data augmentation for some
univariate (r=1) ordinal time series. 
We then employ it to estimate multivariate times series copulas,
where MCMC cannot be used in practice.

Ordinal time series data arise in many fields,
such as criminology <cit.>, marketing <cit.> and finance <cit.>. 
These series often exhibit over- or under-dispersion, multiple modes, truncation and zero-inflation
in the margin, along with serial correlation in the level and also conditional variance (ie. heteroskedasticity). There is an extensive literature
on models that can capture one or more of these
features in univariate series;  see <cit.> and <cit.> for some examples. 
In contrast, extension to multiple series is difficult and less common; for examples,  see <cit.> and <cit.>. 
In comparison, time series copula models
allow for over- or under-dispersion, multiple modes, truncation and zero inflation in a time 
series through the choice
of an arbitrary margin. They also extend readily to multiple ordinal or mixed time series. 
However, a
major challenge is selecting a high-dimensional copula that 
can capture both persistence in the mean and  heteroskedasticity parsimoniously. To do so, we employ a drawable vine (D-vine) <cit.>,
which <cit.> and <cit.> show is parsimonious for Markov and stationary multivariate time series.
Following <cit.>, the component
`pair-copulas' are carefully selected to capture heteroskedasticity in the series.
This is important because it is a feature exhibited by most ordinal time series.

<cit.> use low-dimensional copulas 
to capture cross-sectional dependence between multiple ordinal-valued time series. 
However, this is different from what we propose here, 
where we employ intrinsically high-dimensional copulas to capture both serial and cross-sectional
dependence jointly.  In early work, <cit.> outlined a copula-based Markov time
series model for ordinal data, which <cit.> extend to multivariate panel data using 
a low-dimensional elliptical copula to capture contemporaneous cross-sectional dependence. While this model is
parsimonious and tractable, it does not allow for direct dependence
between lagged values of series. This is often
important in the time series modeling of multivariate continuous data, and we find it is also 
important in our multivariate examples here.







The efficacy of both the VBDA estimator and the proposed copula time series model 
is illustrated using six example time series. The first three are monthly counts of  
murder, attempted murder and manslaughter in the Australian state of New South Wales. 
Univariate time series copula models
show that both counts of murder and attempted murder exhibit
serial dependence, including heteroskedasticity. When compared to the (effectively exact)
posterior computed using MCMC data augmentation, the VBDA estimates prove highly 
accurate, yet are much faster to compute. The fourth example is a 
binary-valued time series 
simulated from an auto-logistic regression with very high serial dependence.
This is an extreme test,
for which MCMC data augmentation fails, yet VBDA gives good results.   
The VBDA estimator is then applied to a trivariate 
time series copula model of the three crime count series.
The 792-dimensional copula  
captures a rich multivariate serial dependence
structure, and is difficult to estimate using MCMC data augmentation in reasonable time.
The last example illustrates the mixed margin case, where 
the bivariate serial dependence structure of monthly counts of U.S. bankruptcies and 
the VIX
(which is a continuous-valued index of stock market volatility) is estimated. The marginal distribution 
of the VIX is highly irregular, making the copula model attractive because it can be
modeled nonparametrically. The estimated 658-dimensional copula  
captures heteroskedasticity in both series, and indicates that the VIX is a leading
indicator of U.S. bankruptcies. These copula models cannot be estimated using the 
VBIL method as outlined in <cit.> in reasonable time.

The paper outline is as follows. Section <ref> outlines copula models for discrete
data, and the D-vine copula for univariate time series. Section <ref> presents
our new VBDA method, including different variational approximations to the augmented posterior,  MCMC data augmentation, and the four univariate examples.
Section <ref> extends the vine copula model to the case of multiple ordinal time series,
and Section <ref> to a mixture of ordinal and continuous series. Section <ref> concludes.




§ COPULA MODEL





 §.§ Copula with Discrete Margins



Following <cit.>, the joint distribution function of a discrete-valued random vector
Y=(Y_1,…,Y_T) can be written as
 
    F( y|θ) = C( u|θ),

where y = (y_1,…,y_T)', u =(u_1,…,u_T)', u_t = G_t(y_t), G_t is the marginal distribution function of Y_t, and C is a T-dimensional copula function that captures all dependence in Y. In the copula modeling literature it is usual 
to select a parametric copula for C, with parameter vector θ.
Because Y_t∈ S for a finite or countably infinite set S, then F at Equation (<ref>) is only uniquely defined
on its sample space <cit.>.
Nevertheless, F
remains well-defined for any given parametric copula function C. Let b_t=G_t(y_t), and a_t = G_t(y_t^-) be the left-hand limit of G_t at y_t, then the corresponding probability mass function is 

    f( y|θ) = Δ_a_1^b_1…Δ_a_T^b_TC(ω|θ) ,

where the difference notation of <cit.> is employed with vector of
differencing variables ω. 
Direct computation of Equation (<ref>) is impractical in higher dimensions because
it involves evaluation of C a total of 2^T times. 
However, following <cit.>, likelihood-based estimation can be undertaken by
introducing a latent vector U=(U_1,…,U_T)', such that (Y',U)' have
augmented density

    f(y,u|θ) = c(u|θ)∏_t=1^T I(a_t≤ u_t< b_t) ,

with copula density c( u|θ) = ∂^T/∂ u_1,…,∂ u_TC( u|θ), and 
the indicator variable I(X)=1 if X is true,
and I(X)=0 otherwise. (Note that an alternative notation to the indicator function
here is δ_y_t(G_t^-(u_t)), 
which is a Dirac mass at G_t^-(u_t), with G_t^- the 
quantile function of Y_t.)
The margin in y of Equation (<ref>) 
is the required mass function at Equation (<ref>). 

When there are multiple independent observations on Y, as with
the cross-sectional and longitudinal datasets considered in <cit.>, 
then the augmented likelihood is the product of Equation (<ref>) over the observations. For
the time series case that is the focus of this
paper, the augmented likelihood is given directly by Equation (<ref>).

Last, we note that throughout the paper we denote copula densities with a `c', and
density/mass functions that are posteriors 
with `p', variational approximations with `q', and all others with `f'.





 §.§ Time Series Copula
 


We consider the case where 
{Y_t} is a strongly stationarity ordinal-valued stochastic process with Markov order p.
Then
G_t is time invariant 
and can be written as G, and
the 
main challenge in using the copula model at Equation (<ref>) is the 
selection of C to capture the serial dependence in the 
series. We note that ordinal time series usually
exhibit persistence in both the mean and variance, so that C should capture this feature. 
To do so, we adopt
a D-vine copula <cit.> with pair-copula components
carefully selected to capture persistence in the first two moments.



In general, a D-vine copula density is equal to
the product of T(T-1)/2 bivariate copula densities called pair-copulas. 
However, when the series has Markov order p, the number of pair-copulas is
much smaller. Moreover, when the series is also stationary, the number of 
unique pair-copulas is equal to the Markov order p <cit.>.
For s<t, by denoting u_t|s = F(u_t|u_s,…,u_t-1),
u_s|t = F(u_s|u_s+1,…,u_t)
and u_t|t = u_t,
this parsimonious D-vine copula density is

    c^DV(u|θ)    =    ∏_t=2^T f(u_t|u_max(1,t-p),…,u_t-1) 
       =    ∏_t=2^T∏_k = 1^min(t-1,p)
    c_k+1(u_t-k|t-1,u_t|t-k+1;θ_k+1) ,

where θ = {θ_2,…,θ_p+1} and c_2,…,c_p+1
are the pair-copula densities. Given u,
 the arguments {u_t|s,u_s|t ; t=2,…,T , s<t} 
are computed using the recursive Algorithm 1 in <cit.>.

<cit.> show that c^DV is able to capture persistence
in the variance if one or more c_k allows for concentration of the probability mass in the four quadrants of the unit square. To do so they suggest
 the following mixture of rotated copulas:

    c^MIX(u,v;γ) = wc^a(u,v;γ^a)+(1-w)c^b(1-u,v;γ^b)
     , 0≤ w ≤1 .

Here, γ={γ^a,γ^b,w}, 
0≤ w ≤ 1 is a weight, and c^a,c^b are two parametric bivariate copula densities with non-negative Kendall's tau and parameters γ^a and γ^b respectively. 
In our empirical work, for the mixture components c^a and c^b we employ the 
`convex Gumbel' defined as follows.
Let c^G(u,v;τ) be the density of a Gumbel copula parameterized (uniquely)
in terms of its Kendall tau value 0≤τ < 0.99. (Note that we bound
τ away from 1 to enhance numerical stability of the D-vine copula.) 
Then the convex Gumbel has a density c^cG equal to the convex
combination of that of the Gumbel and its rotation 180 degrees (ie. the
survival copula), so that

    c^cG(u,v;τ,δ)=δ c^G(u,v;τ) + (1-δ)c^G(1-u,1-v;τ) ,

with 0≤δ≤ 1.
When employed for c^a and c^b in Equation (<ref>), it gives a
five parameter bivariate copula with γ^a=(δ^a,τ^a), 
γ^b=(δ^b,τ^b),
and a density c^MIX that is equal 
to a mixture of all four 90 degree rotations of the Gumbel copula. 
We use independent uniform priors on the elements of γ 
in our empirical work. 

To measure 
the level of serial dependence captured by our copula model, we use the Spearman's correlation
between Y_s and Y_t for s<t. Following
<cit.>, for ordinal-valued variables this is

    ρ_k     =3∑_y_s∈ S∑_y_t∈ Sg(y_s)g(y_t) (C̅_k(b_s,b_t) +C̅_k(b_s,a_t)  + C̅_k(a_s,b_t) +C̅_k(a_s,a_t)) -3 ,

where C̅_k is the copula function of the distribution of 
(Y_s,Y_t), which only varies with k=t-s when Y_t is stationary <cit.>. 
This copula C̅_k is constructed by 
simulating (many) draws of u from c^DV using Algorithm 2 in <cit.>, and then constructing the bivariate empirical
copula from the draws of elements (u_s,u_s+k).







§ BAYESIAN ESTIMATION



From Equation (<ref>), the augmented posterior density is

    p(θ,u|y) = f(y,u,θ)/f(y)=( c(u|θ)p(θ)∏_t=1^T I(a_t≤ u_t< b_t) )/f(y) ,

where p(θ) is the prior and f(y) is the
 marginal likelihood. The augmented posterior above admits p(θ| y) as one of its margins, and is 
tractable (up to proportionality). <cit.> propose a MCMC data augmentation method for its (effectively exact) evaluation. However, 
this MCMC scheme is generally slow, and computationally infeasible for 
high-dimensional copulas.
Variational Bayes (VB) is an alternative inferential method to MCMC, with both methods typically applicable to the same problems.
Here, we use the augmented posterior above to develop a new
VB estimator for p(θ| y). 





 §.§ Variational Bayes Estimator



VB makes possible the estimation of copula models with discrete margins, even for copulas in high dimensions and with a large number of parameters.
Here, p(θ,u|y) is approximated by a tractable density  q_λ(θ,u) with parameters λ, called the variational
approximation. 
Estimation consists of finding values of λ that minimize the Kullback-Leibler divergence  

    KL(q_λ(θ,u)||p(θ,u|y)) =∫log(q_λ(θ,u)/p(θ,u|y))q_λ(θ,u)θu .

This can be shown <cit.> to correspond to maximizing the lower bound of the logarithm of the marginal likelihood log p( y), given by

    ℒ(λ)=∫log(p(θ)
    	f(y,u|θ)/q_λ(θ,u))q_λ(θ,u)θu .

In selecting q_λ, it is common to assume independence between some or all
parameters  <cit.>, and we do so here 
between θ and U. The variational approximation we use has density

    q_λ(θ,u)=
    q_λ^a(θ) q_λ^b(u)
     ,

where the density
q_λ^a has parameters λ^a,
the density q_λ^b has support on [a_1,b_1)×…× [a_T,b_T) and parameters λ^b, and λ={λ^a,λ^b}.
The key to the success of our method is the specification of 
q_λ^a and q_λ^b, which we discuss in detail later.

We follow <cit.> and
others and use
stochastic gradient ascent (SGA) methods to maximize
 ℒ(λ).
 This approach only requires that (i) generation from q_λ(θ, u) 
 is possible, and that (ii) the target distribution is tractable and can be evaluated up to proportionality. 
 Condition (i) is met by our choices for  q_λ^a and q_λ^b outlined below.
 Condition (ii) is met because the augmented posterior is 
 tractable, whereas p(θ|y) based on Equation (<ref>) is
 not. To implement SGA, initial values for the parameters, λ^(0), are selected and then the lower bound is sequentially optimized by
 values λ^(1),λ^(2),… obtained by the updating formula

    λ^(k+1) = λ^(k)+ρ^(k)∇_λℒ(λ^(k)) .

Here, ∇_λℒ(λ^(k)) is an unbiased estimate of the lower bound's gradient 
 ∇_λℒ(λ^(k)), and ρ^(k) is the learning rate, set using the ADADELTA method described in the Appendix <ref>.
To compute the gradient, SGA methods resort to the “log-derivative trick” (E_q(∇_λlog q_λ(θ, u)) = 0), and show that 

    ∇_λℒ(λ)=E_q( ∇_λlog q_λ(θ,u){log h(θ,u)-log q_λ(θ,u)}) ,

with p(θ,u|y)∝ f(y,u|θ)
p(θ)=h(θ,u), and E_q is the expectation with respect to q_λ(θ,u). Notice from Equation (<ref>) that an unbiased estimate is  
∇_λℒ(λ)= (g_λ_1,…,g_λ_m)', where

    g_λ_i=1/S∑_s=1^S(log h(θ_s,u_s)-log q_λ(θ_s,u_s))∇_λ_ilog q_λ(θ_s,u_s) ,

with m as the number of elements in λ. 
An advantage of our choice of variational approximation is that 
the ith element of the gradient ∇_λlog q_λ(θ,u)
simplifies to
∇_λ_ilog q_λ(θ,u)= 
∇_λ_ilog q_λ^a(θ)+ 
∇_λ_ilog q_λ^b(u).
For s=1,…,S, the values 
θ_s∼ q_λ^a(θ) and u_s=(u_1,s,…,u_T,s)∼ q_λ^b(u) 
on [a_1,b_1)×⋯× [a_T,b_T). 

Algorithm <ref> presents how the SGA optimization works within variational Bayes. Step (1b) is based on the work by <cit.>, which employs a vector of control variates, ς, for variance reduction of the unbiased estimate of the gradient. The stopping rule is commonly set as a fixed number of SGA steps taken <cit.>.






 §.§ Variational Approximation



Key to developing an effective VB estimator is the selection of q_λ^a and q_λ^b in 
Equation (<ref>) that balance tractability and accuracy. 
We first outline
three choices for q_λ^b, after which we then detail that for q_λ^a.





  §.§.§ Approximation for U



To guide our choice for q_λ^b, we derive the marginal posterior of
U in Theorem <ref>.

If (θ,U) have the augmented posterior density function at Equation (<ref>), then: 

	
  (a) The joint density p(u|y)=c̃(u) ∏_t=1^T  I(a_t≤ u_t < b_t)/f(y),
	where c̃(u)=∫ c(u|θ)p(θ)θ is a copula density and f(y) is the marginal likelihood; and,
	
  (b)
	the marginal density p(u_t|y)∝ I(a_t≤ u_t < b_t)∫ A(u_t|θ)p(θ|y)θ, where 
	
    A(u_t|θ)=Δ_a_1^b_1⋯Δ_a_t-1^b_t-1Δ_a_t+1^b_t+1⋯Δ_a_T^b_T H(v_1,…,v_t-1,u_t,v_t+1,…,v_T|θ) ,

	H(u|θ)=∫ c(u|θ)u_s≠ t  and u_s≠ t=(u_1,…,u_t-1,u_t+1,…,u_T).

Proof: See Appendix <ref>.

We make two observations on the  posterior of U.
First, if the elements of Y are independent, then c̃(u)=1 is the density 
of an independence copula, and p(u|y)∝∏_t=1^T  I(a_t≤ u_t < b_t), so that each element is independent uniform. Second, as (b_t-a_t) → 0 for 
all t, then p(u|y)∏_t=1^T δ_u_t(a_t). That is, as the data Y becomes
`closer to continuous', the posterior approaches a degenerate distribution with point mass at  u=(a_1,…,a_T)'.
 
Armed with these observations, our first choice is simply independent uniforms:

    q_λ^b(u)=∏_t=1^T 1/b_t-a_t I(a_t ≤ u_t < b_t) , λ^b=∅ .

We expect VA1 to be more accurate for data with low dependence (although we find it still 
works well for even quite dependent data).
The next two approximations are based on normal distributions for a transformation 
of U. Let
Z_t=Φ^-1((U_t-a_t)/(b_t-a_t)), and
Z=(Z_1,…,Z_T)'∼ N(η,Ω), with
η=(η_1,…,η_T)' and Φ the standard normal distribution
function. The Jacobian of this transformation is 
J_Z→U=∏_t=1^T ((b_t-a_t)ϕ(z_t))^-1, with
ϕ
the standard normal density and z_t=Φ^-1((u_t-a_t)/(b_t-a_t)).
Our second choice for q_λ^b assumes 
Ω=(ω^2_1,…,ω^2_T), so that

    q_λ^b(u)=
    	∏_t=1^T ϕ_1(z_t;η_t,ω^2_t)/(b_t-a_t)ϕ(z_t) ,λ^b={η,logω_1,…,logω_T} ,

and ϕ_1(z_t;η_t,ω^2_t) is the density of a N(η_t,ω^2_t) distribution
evaluated at z_t. Note that VA2 nests VA1.
We find this an effective mean field approximation that is accurate for a wide range of data,
and very fast to work with. 

For Y that exhibits extreme dependence,
our third choice allows for the elements
of U to be dependent[We 
are grateful to an anonymous
referee who suggested that this may be an important consideration.]
by adopting a non-diagonal (but sparse) precision matrix Ω^-1.
For time series copulas, we set Ω^-1=LL', with L a band one lower triangular
Cholesky factor. This corresponds to an approximation q_λ^b with the dependence
structure of a (non-stationary) first order Markov process for {U_t}. The density

    q_λ^b(u)=
     ϕ_T(z;η,(LL')^-1)/∏_t=1^T(b_t-a_t)ϕ(z_t) ,λ^b={η,L} ,

where only the free elements of L are variational parameters. If the lower triangular first band
of L contains only zeros, then VA3 reduces to VA2. Both approximations are fast to generate
from by first generating z from normals, and then transforming to 
u. The gradients ∇_λ^blog q_λ^b(u) required to 
implement Steps 1(b) and 2(b) of Algorithm <ref> are available in
closed form; see Appendix <ref>.

Note, as (b_t-a_t) → 0 for all t, all three approximations 
become exact. The accuracy of q_λ^b is important because
it can also increase the accuracy of the variational approximation of p(θ|y).





  §.§.§ Approximation for θ



Denoting the number of parameters in θ as n, the most popular choice for 
q_λ^a is the density of a N(μ,Σ) distribution, because it is quick to
generate from and the gradient ∇_λ^alog q_λ^a(θ)
is available in closed form ().

To ensure Σ is positive definite, λ^a is typically a convenient
re-parametrization of μ and Σ. 
In applications where θ has a large number of elements, a sparse representation of Σ helps to improve the accuracy of the gradient estimate and its speed of computation. We follow <cit.>, who suggest the factor representation of the covariance matrix Σ = B'B+D, where the matrix B is of dimension n× K, K is the number of factors and K<<n. All the elements in the upper triangle of B are set to zero. D is a n× n diagonal matrix such that D_i,i = d_i^2, where d_i is the i^th element of the vector d. <cit.> derives the gradient for this case, and shows it is
fast to compute; see also Appendix <ref>. In our empirical work, we
compare the accuracy of the 
approximations for various values of K and find low values adequate.





  §.§.§ Discussion of Alternative VB Approximations



<cit.> suggest using an unbiased estimator of the
intractable likelihood in Equation (<ref>) computed using importance sampling. This involves
drawing N_IS values of u, at which c(u|θ) is repeatedly evaluated.
Whenever evaluating the copula density is computationally intensive — such as for c^DV here or with 
other high-dimensional or complex copulas —  this will be many times slower than our approach.

Variational approximations to posteriors augmented with latent variables
have proven successful in a number of other models; see <cit.>, <cit.> and <cit.>
for some recent examples. However, <cit.> and others observe that assuming independence 
between the latent variables in such an approximation may lead to poor inference in some circumstances. This motivates
VA3, although in our
empirical work VA2 proves almost as accurate and several times faster.
Key to using VA3 for other copula models is the adoption of
an appropriate parsimonious matrix Ω, or its inverse. Last, we mention
it is also possible to employ
a Gaussian approximation with factor covariance 
structure for the vector (θ,Z).  This may improve the accuracy of the approximation
for some copulas, but will introduce an extra KT variational parameters
(ie. the extra factor loadings),
slowing estimation down substantially. Our empirical work suggests that our proposed variational approximations strike a balance between computation
time and accuracy.





 §.§ Data Augmentation



We now outline MCMC data augmentation, tailored for the parsimonious
D-vine copula in Section <ref>. Key to implementation is
the evaluation of the conditional
densities and distribution functions below.
If t_0=max(t-p,1) and 
t>1, then

    f(u_t|u_t_0,…,u_t-1)    =    ∏_k=1^min(t-1,p)c_k+1(u_t-k|t-1,u_t|t-k+1;θ_k+1) 
    
    	F(u_t|u_t_0,…,u_t-1)    =    h_t_0,t∘ h_t_0+1,t∘⋯∘ h_t-1,t(u_t) ,
  
where h_s,t(u)=h_t-s+1^1(u|u_s|t-1)=∂/∂ vC_t-s+1(v,u)|_v = u_s|t-1 is the conditional pair-copula function. This is given in
Appendix C1 of <cit.> for the mixture copula
defined at Equation (<ref>).

The values u are integrated out of the augmented posterior
as part of an MCMC sampling scheme. The scheme generates from the conditional
posteriors (1) p(u|y,θ),
and (2) p(θ|u). 
Given the values u, step (2) can be undertaken using (adaptive) random walk Metropolis-Hastings (MH), where
θ_k is generated conditional on 
{θ\θ_k} for k=2,…,p+1.
However, step (1) is more involved, with the latent variables u 
generated jointly using a MH
step. <cit.> suggest using the
proposal density
π(u) = ∏_t=2^Tπ_t(u_t|u_t_0,…,u_t-1)π_1(u_1),
where
π_1(u_1) = ℐ(a_1≤ u_1<b_1)/(b_1-a_1)
and

    π_t(u_t|u_t_0,…,u_t-1) ∝f(u_t|u_t_0,…,u_t-1)ℐ(a_t≤ u_t<b_t)/F(b_t|u_t_0,…,u_t-1)-F(a_t|u_t_0,…,u_t-1) .

Therefore,
a proposal iterate
u^ can be
obtained from π(u) by generating sequentially from
the univariate densities π_1,…,π_T. Each of these is a
constrained univariate distribution with known distribution function, so that
iterates can be generated easily using the inverse distribution 
method. An advantage of the proposal π is that
the
MH acceptance ratio is fast to compute. The probability of
accepting 
u^=(u_1^,…,u_T^) over the previous value 
u^=(u_1^,…,u_T^) is

    min(1 , ∏_t=2^TF(b_t|u_t_0^,…,
    	u_t-1^)-F(a_t|u_t_0^,
    	…,u_t-1^)/F(b_t|u_t_0^,…,u_t-1^)-
    	F(a_t|u_t_0^,…,u_t-1^)) .


In general, this proposal works well. However, for challenging high-dimensional
copulas with highly dependent binary-valued data Y, we find that the acceptance
rate for this step can be prohibitively low and the MCMC scheme can get stuck. We illustrate
this empirically below. Throughout, we
employ a burnin sample of 10,000 iterates, followed by 
a further 20,000 iterates from which we compute posterior inference.

 
 
	              





 §.§ Examples



We illustrate the copula time series model, and the efficacy of the VB estimator, using four univariate time series examples. The first three are T=264 monthly counts of the crimes of Murder, Attempted Murder and Manslaughter
in the Australian state of New South Wales (NSW) between January 1995 and December 2016. The data is sourced from the NSW Bureau of Crime Statistics and Research. The fourth is T=200 observations simulated from an 
auto-logistic regression with
(Y_t=1|y_t-1)=^-1(-2.197+4.394y_t-1).
This last example is challenging because (i) the series is short, (ii) the data are binary-valued, (iii) there is very strong serial dependence, and (iv) there are
many more 1's. 
Figure <ref> gives the relative frequency histograms and time series plots of the four series.




We set p=3, and fit the copula using c^MIX pair-copula components for c_2, c_3 and c_4. For c^a and c^b we chose convex Gumbels, so that θ_k = {τ_k^a,δ_k^a,τ_k^b,δ_k^b,w_k} for k = 2,3,4. This T-dimensional
D-vine copula has a total of n = 15 parameters.
We set G to the empirical distribution functions in Figure <ref>. 
To estimate the copula parameters we first use MCMC data augmentation to compute the exact
posterior as outlined in Section <ref>.
For the three crime series, between 22% and 71% of MH proposals were accepted, but
for the Auto-Logistic example
the sampler became stuck, and estimation failed.
For the crime series,
Table <ref> reports the posterior means and intervals
 of the copula parameters. To summarize the serial dependence captured by the copula, Table <ref> reports the posterior of the Spearman correlations
 ρ_k for k=1,…,3
 for the three crimes. Correlation is strong for Attempted Murder, but not for Murder and 
 Manslaughter. However, this measures correlation in the level of the series, and not more general dependence, such as in higher order moments. Figure <ref> presents the log-densities of the pair-copulas at the posterior mean values. Most of these copula densities are far from uniform, indicating more general serial dependence exists in these series. The pair-copulas have 
 probability mass in the off-diagonal corners of the unit square, which <cit.>
 show is indicative of serial correlation in conditional variance (ie. heteroskedasticity).
 
 
 
 
  
 

The four D-vines were also estimated using VBDA. 
Because all parameters are bounded between 0 and 1,
we transform them to the real line as θ_k =
 {ψ(τ_k^a),ψ(δ_k^a),ψ(τ_k^b),ψ(δ_k^b),ψ(w_k)},
 where ψ(a) = log(a/1-a). Estimation was implemented separately 
for approximations VA1, VA2 and VA3, with K = 0,…,15 factors. Each
estimator used 5000 SGA steps, and S=500 to 
estimate the gradient.
The initial values for λ are
{B=0, D=√(0.1)I_n, μ=μ_0, η=0, L=I_n}, where
μ_0 is set to values where the pair-copulas are all independence pair-copulas. 
These initial values are used in all our empirical work, although the results are robust to changes in them.
Figure <ref>(a,c,e,g) shows how the lower bound increases with K, and any increase is small for K≥ 3.
Figure <ref>(b,c,f,h) plots the lower bound against SGA step when K=3,  suggesting that the SGA algorithm converges within 1000 steps in every case. 
For the three crime series, the lower bounds of 
VA2 and VA3 are almost indistinguishable. For the challenging Auto-Logistic example, 
VA3 — which allows for dependence in the latent variables — dominates. However, it is 
difficult to determine how much of the higher lower bound values are attributable 
to an increase in the accuracy of q_λ^a, as opposed to q_λ^b. 
Nevertheless, plots of the pair-copula densities (see Supplementary Materials) for each of
VA1, VA2 and VA3 suggest that VA3 provides a meaningful 
improvement over VA1 and VA2. These also show the VBDA estimates suggest the series
has Markov order one, and high serial dependence; which correspond to the known data generating
process.
 
 

To illustrate the accuracy of VBDA for the three crime series, 
Figure <ref> plots the posterior means and standard deviations of θ from the preferred approximation (VA2 with K=3)
against their (effectively exact) values computed via MCMC. 
Both moments of the VB approximations are close to those of the true posterior.
Similar plots for VA1 and VA3 (see Supplementary Materials)
suggest these are also reasonable approximations. The
first three rows in Table <ref> present the copula specifications and total estimation times for MCMC and VBDA for all examples. The computations
were undertaken on a Dell Precision workstation using Matlab, and in parallel using
8 workers for key computations for both estimators. The results show that VBDA is many times faster than
MCMC data augmentation. Moreover, the main computation of the VBDA estimator
is the repeated evaluation of  
h at Step 2(b). This is slow because computing the arguments of the pair-copulas
is computationally intensive, and the
VBDA method proves even faster for 
simpler copulas.



 





§ MULTIVARIATE ORDINAL TIME SERIES



In this section we extend the time series copula to capture the dependence in 
multiple ordinal-valued series.





 §.§ Copula Model and Estimation



Consider an r-dimensional stationary stochastic process {Y_t}, where 
Y_t=(Y_1,t,…,Y_r,t)', and each element Y_i,t is ordinal-valued
with margin G_i. We assume
a time series copula model with Tr-dimensional copula function.
Then if y = (y_1',…,y_T')'
and 
y_t = (y_1,t,…,y_r,t)', we estimate it using the
augmented likelihood

    f(u,y|θ) = c(u|θ)∏_i=1^r∏_t=1^T I(a_i,t≤ u_i,t< b_i,t) ,

where u = (u_1',…,u_T')',
u_t = (u_1,t,…,u_r,t)', a_i,t=G_i(y_i,t^-) and b_i,t=G_i(y_i,t).


The copula density c in Equation (<ref>)
 captures both cross-sectional and serial dependence jointly. For this, 
<cit.> and <cit.> use a Gaussian copula,
with parameter matrix equal to the correlation matrix of a stationary vector autoregression. However, a Gaussian copula cannot capture the high level of persistence in the variance often
exhibited in ordinal time series. Instead, we follow <cit.>, <cit.>, <cit.> and <cit.> and again use a D-vine copula, but with a parsimonious
form corresponding to a stationary Markov p multivariate series. 
The pair-copula components are of the form at Equation (<ref>) to account for heteroskedasticity.

<cit.> shows that this D-vine has a density that can be factorized as

    c^DV(u)= K_0(u_1)∏_t=2^T(  K_0(u_t) 
    ∏_k=1^min(t-1,p) K_k(u_t-k,…,u_t)) .

The functionals K_0,…, K_p are each products of blocks of
pair-copula densities,
and do not vary with t for stationary series. They are
defined as

    K_k(u_t-k,…,u_t) = 
    ∏_l_1=1^r∏_l_2=1^l_1-1c_l_2,l_1^(0)(u_j|i-1,u_i|j+1;θ_l_2,l_1^(0))   if  k=0
    ∏_l_1=1^r∏_l_2=1^rc_l_2,l_1^(k)(u_j|i-1,u_i|j+1;
    θ_l_2,l_1^(k))   if  1≤ k≤ p ,

where c^(k)_l_2,l_1 is a bivariate pair-copula density 
with parameters θ_l_2,l_1^(k). 
When k=0, there are r(r-1)/2 of these associated with K_0, 
and they collectively capture
cross-sectional dependence between the r variables. For example,
if they were each equal to the bivariate independence copula with 
density c^(0)_l_2,l_1=1, then 
K_0=1 and the variables would be independent
contemporaneously. When k>p, there
are r^2
pair-copulas associated with block K_k that capture serial dependence 
at lag k. 
In total,
there are p(r^2)+r(r-1)/2 unique pair-copulas, which is much less than
the Tr(Tr-1)/2 in an
unconstrained D-vine. 
The indices of the pair-copula 
arguments are i=l_1+r(t-1) and j=l_2+r(t-k-1), 
and
the argument values
{u_i|j,u_j|i; i=1,…,Tm, j<i} 
are computed
using the Algorithm 1 of <cit.>.
Last, we note that if r=1, then K_0=1, i=t, j=t-k and 
K_k=c_1,1^(k)(u_t-k|t-1,u_t|t-k+1), so that with the notation
c_k+1≡ c_1,1^(k), the copula densities at 
Equations (<ref>) and (<ref>) are the same.

To measure the dependence between Y_j,s∈ S_j and Y_i,t∈ S_i, with k=t-s, we use the Spearman's correlation

    ρ_i,j,k =-3+3∑_y_j,s∈ S_j∑_y_i,t∈ S_i   g_j(y_j,s)g_i(y_i,t)(C̅_j,i,k(b_j,s,b_i,t).+C̅_j,i,k(b_j,s,a_i,t)  +
       .  C̅_j,i,k(a_j,s,b_i,t) +C̅_j,i,k(a_j,s,a_i,t)) .

Here, g_i is the probability mass function corresponding to G_i,
while C̅_j,i,k(u_j,s,u_i,t) is the copula function
of the bivariate marginal of (Y_j,s,Y_i,t). The latter is computed by simulating from
 c^DV and then constructing the empirical copula function for 
C̅_j,i,k.

The augmented posterior of this copula time series model is

    p(u,θ|y) ∝ c^DV(u|θ)p(θ)∏_i=1^r∏_t=1^T I(a_i,t≤ u_i,t< b_i,t) .

Because of the very large number of elements in u, 
estimation using MCMC is computationally infeasible for  
even moderate values of r and T.
However, our VBDA estimator can be employed with the same variational approximations
outlined in Section <ref>. We note that in our empirical work we employ VA3 as exactly outlined, although
the sparse pattern of
Ω^-1 can be further tailored to match the possible dependence
structure of p(u|y) for this case.








 §.§ Example: New South Wales Homicide



We consider a trivariate time series copula model for the NSW monthly crime counts, with 
the empirical distributions as univariate marginals.  
The copula density is given in
Equation (<ref>), where we set p=1 and adopt 
pair-copula densities of the form c^MIX. 
The dimension of the D-vine copula is 3 × 264 = 792, 
and Table <ref> reports its specification. The copula parameters are 
estimated using the VB estimator with
K = 0,1,2,4,5,10,15,40,50 factors.  
To estimate the gradient S=500 in Algorithm 1, and 5000 VB steps are used with K=15.
Figure <ref> plots the variational lower bound 
against K in panel (a), and against the VB step when K=15 in panel (b), for VA1, VA2 and VA3.
 A total of 
K=15 factors appears sufficient, while both VA2 and VA3 give similar results, but dominate VA1. 
Table <ref> reports the posterior means 
and intervals of θ for VA2, although those for VA1 and VA3 are very similar and are reported
in the Supplementary Material.

      



Table <ref> reports the
estimates of the pairwise Spearman correlations. The 
contemporaneous correlations (k=0) are given in the top two rows, and first order serial correlations
(k=1) in the bottom rows.
There is positive 
contemporaneous correlation between Attempted Murder and Murder, and also (weakly) with Manslaughter. There is first order serial correlation in
Attempted Murder, but not in the other two crimes. The most striking result is that 
Attempted Murder is positively correlated with Murder and Manslaughter one month later, suggesting
it is a leading indicator of these two crimes.
However, these correlations measure dependence in the level only.  
Figure <ref> displays the logarithm of
the 15 unique pair-copula densities. 
Most have mass in the off-diagonal corners of the unit square, indicating
that the copula is capturing heteroskedasticity and `variance spill-overs'
between the three series. The 3 pair-copulas on the lefthand side capture
contemporaneous cross-sectional dependence. The 9 pair-copulas on the righthand side 
capture first order serial dependence. For example, c_2,1^(1) is very far from
uniform, and captures strong variance spill-over between Attempted Murder and Murder.







§ MIXED MULTIVARIATE TIME SERIES





 §.§ Copula Model and Estimation



Consider the case of a stochastic process {Y_t}, where
Y_t=(Y_1,t,…,Y_r,t)' consists
of d ordinal and r-d continuous-valued variables, which we refer to as `mixed'. A copula model using the same D-vine with density 
c^DV at Equation (<ref>) can be used 
to construct a flexible time series model. 
Without loss of generality,
if the first d elements of Y_t are ordinal,



then <cit.> suggest estimation of a copula model with mixed margins using the 
augmented density at Equation (<ref>), but
where f(y_i,t|u_i,t)= I(u_i,t=G_i(y_i,t)) is a point
mass for i=d+1,…,r. They
discuss how to implement MCMC data augmentation,
but this approach can be slow or computationally
infeasible for values of Td that occur frequently in time series analysis.

Let C={(i,t): i=d+1,…,r ; t=1,…,T} denote
the indicies of the continuous-valued Y_i,t, and 
u_D be the Td latents corresponding
to the ordinal variables. Then 
VBDA can be employed using the variational
approximations outlined in Section <ref>, but where
q_λ^b(u)=
q̃_λ^b(u_D)∏_(i,t)∈ C I(u_i,t=G_i(y_i,t)), and 
approximations VA1 to VA3 are considered for q̃_λ^b. 
Algorithm <ref> can be used to approximate
the augmented posterior, but where u_i,t=G_i(y_i,t) are constants
for (i,t)∈ C, and are not generated.

Equation (<ref>) can be used to compute the Spearman correlation ρ_i,j,k  between two 
ordinal-valued variables (Y_j,s,Y_i,t) with k=t-s and s<t. If  both variables are 
continuous-valued, then
ρ_i,j,k=12∫C̅_i,j,k(u,v)uv-3. But if Y_j,s∈ S_j is ordinal and Y_i,t is
continuous, then 

    ρ_i,j,k=
    6∑_y_j,s∈ S_jg_j(y_j,s)∫ g_i(y_i,t)(C̅_j,i,k(b_j,s,G_i(y_i,t)) + C̅_j,i,k(a_j,s,G_i(y_i,t)))y_i,t -3 ,

where the integral can be computed numerically. In all cases, C̅_i,j,k
is evaluated by simulation as previously.



          







 §.§ Example: Bankruptcy and the VIX



We study the dependence between the continuous-valued
VIX index, which measures
U.S. market volatility, and the number of public company bankruptcy cases filed
in U.S. courts.
We employ the monthly average value of the VIX
obtained from the FRED website, while the bankruptcies are monthly counts
sourced from the UCLA-LoPucki Bankruptcy Research Database. The time series
are 
from
December 1989 to April 2017, so that T=329.

A positive relationship between market volatility and bankruptcies has been documented previously   <cit.>. Figure <ref>(a) plots both series, 
while Figure <ref>(b) displays the empirical
distribution of the VIX, conditional on the number of bankruptcies.
The positive correlation between the two series
is apparent in both panels.

We employ the D-Vine copula model with p=2 and a total number of 5×9 = 45 parameters;
Table <ref> reports the copula
specification. Separate 
variational approximations with K = 0,1,2,4,5,10,15,40,45 factor decompositions
for Σ were estimated. The 
same values for S and number of VB steps were adopted as in Section <ref>.
Figure <ref>(a) plots the lower bound against K, and it varies
little for K≥ 15. This is consistent with the empirical results in <cit.>, who
found that a higher number of factors are needed for the accurate approximation of more complex posteriors.

The copula parameter estimates are reported in the Online Appendix, while
Table <ref> reports the 
pairwise Spearman correlations. Both the number of bankruptcies and the VIX exhibit
serial correlation, although the latter more so. The two series are positively correlated, both 
contemporaneously and in the lagged values.
However, the lagged values of the VIX are more highly correlated with later bankruptcies,
suggesting that the VIX is a leading indicator of public company bankruptcy filings. 
The fitted pair-copulas densities are plotted
in the Online Appendix, and their form
is consistent with heteroskedastic time series. For example, c^(1)_2,1 and 
c^(1)_1,2 have mass concentrated in all four corners, indicating
positive cross-correlation in the variance of the two series at different lags;
ie. volatility `spillover'. 





§ DISCUSSION



This paper makes two main contributions. The first is to propose a new
VB estimator for copula models with discrete, or a combination of discrete
and continuous, variables. The approach can be used to estimate
copulas with a higher dimension and number of parameters than previous methods. We illustrate this using time series copulas of
up to 792 dimensions and 60 parameters, although the method can be used to estimate
copula models for cross-sectional, longitudinal or spatial data 
just as readily.
The second main contribution of the paper is
to propose a new time series model for multivariate ordinal-valued variables, where
a copula captures serial and cross-sectional dependence jointly. 
Our proposed
copula is a parsimonious
D-vine that can capture serial dependence in both the level and 
conditional variance, with the latter being an important feature
in much ordinal data. The time series
model is highly flexible, where any marginal features in a time
series can
be captured by an arbitrary distribution, and is easily extended to a combination of 
discrete and continuous-valued series.
 
Ordinal time series frequently
exhibit both serial dependence in the level and heteroskedasticity.
Few existing copulas can capture both jointly, yet the D-vine used here
can do so for Markov and stationary series. 




An advantage of such copula
time series models is that they allow for the more accurate modeling of data with 
multi-modal and irregular margins, 
as well as being readily extended to multivariate series. However,
their estimation is computationally challenging using previous methods, 
and our VBDA estimator provides a new and effective solution, as illustrated by our
examples.

In our VB approach, a key observation is that it is computationally advantageous to employ
a variational approximation to the augmented posterior p(u,θ|y),
rather than the intractable posterior p(θ|y). This is consistent with
<cit.>, <cit.> and <cit.>, who also find that 
variational approximations to the posterior of latent variables can also 
be computationally efficient in mixed effects 
generalized linear models.
The empirical examples illustrate that the
 approximations proposed here provide a balance between computational
efficiency and accuracy.

Last, we outline some promising directions for future research.
First, <cit.> give an augmented likelihood
for copula models where the margins have mixed densities (not to be confused
with a combination of continuous and discrete variables). 
Extending our VBDA approach to this case would provide a
faster estimator than MCMC.
Second, copula models
for discrete spatial data are growing in popularity
<cit.>. However, estimation is challenging 
for a large number of spatial locations, and
VBDA provides a solution. When employing 
VA3, Ω can be tailored to each case; for example,  
a natural choice for Ω^-1 for data located 
on a regular lattice is the precision matrix of a Gaussian Markov random field.
Last,  copulas constructed by inversion of existing distributions
are popular, including those that have intractable copula 
functions C and densities c; see <cit.> and <cit.> for examples. Extending our VBDA estimator to such intractable copula
models for 
discrete data is an interesting extension.





§ ADADELTA



The learning rate ρ^(k) can be set using different methods. For example <cit.> set it as a sequence with manually tuned parameters. <cit.> propose an adaptive learning rate based on previous work by <cit.>. Here, we employ the ADADELTA method of <cit.>,  which provides reliable convergence of the SGA algorithm. 
This method consists of individually updating the step size for each element in λ as

    λ_i^(k+1) = λ_i^(k)+Δλ_i^(k)

with Δλ_i^(k) = ρ_i^(k)g_λ_i^(k), g_λ_i^(k) denoting the i^th element of ∇_λℒ(λ^(k)) and ρ_i^(k) is given by

    ρ_i^(k)=√(E(Δ_λ i^2)^(k-1)+ϵ)/√(E(g_λ_i^2)^(k)+ϵ)

where ϵ is a small scalar and E(Δ_λ i^2)^(k) and E(g_λ_i^2)^(k) are recursively updated as

    E(Δ_λ i^2)^(k)=   ζ E(Δ_λ i^2)^(k-1)  +(1-ζ)Δλ_i^(k)^2
    
    E(g_λ_i^2)^(k)=   ζ E(g_λ_i^2)^(k-1) +(1-ζ)g_λ_i^(k)^2

For the VB applications here, <cit.> is followed, and we set ϵ =10^-6, ζ=0.95, E(Δ_λ i^2)^(0)=0 and E(g_λ_i^2)^(0)=0.





§ PROOF OF THEOREM 1



The expression at (a) is obtained simply by integrating over Equation (<ref>)
with respect to θ. To derive the expression at (b), first note that 
p(u|θ,y)∝∏_t=1^T  I(a_t≤ u_t < b_t)c(u|θ), so that

    p(u_t|θ,y)    =   ∫ p(u|θ,y)u_s≠ t
       ∝    I(a_t≤ u_t < b_t) ∫ c(u|θ) ∏_s≠ t I(a_s≤ u_s < b_s)u_s≠ t
       =     I(a_t≤ u_t < b_t) A(u_t|θ) ,

where A(u_t|θ) is as defined in Theorem 1.
Therefore, 

    p(u_t|y)=∫ p(u_t|θ,y)p(θ|y)θ∝ I(a_t≤ u_t < b_t) 
    ∫ A(u_t|θ)p(θ|y)θ ,

which is intractable. 
 





§ DERIVATIVES



In this appendix we compute the gradient
∇_λlog q_λ(θ,u) to implement Steps 1(b) and 2(b) for Algorithm <ref>.
To  present these succinctly the following notation is introduced. For a matrix A of dimension n× K, the function vech(.) is defined as vech(A)= (A_1:n,1',…,A_K:n,K')' with A_k:n,k = (A_k,k,…,A_n,k)' for k = 1,…,K. Also, the vector of diagonal entries of the square matrix Z is written as diag(Z). 
Employing this notation, the vector of parameters λ^a can be written as λ^a = (μ',b',d')' with b = vech(B), and the gradient 
∇_λ^alog(q_λ^a(θ))=
(∇_μlog(q_λ^a(θ))',
∇_blog(q_λ^a(θ))',
∇_dlog(q_λ^a(θ))')'
where

    ∇_μlog(q_λ^a(θ)) =   (B'B+D)^-1(θ-μ)                    
    ∇_blog(q_λ^a(θ)) =   vech(-(B'B+D^2)^-1B+(B'B+D^2)^-1(θ-μ)(θ-μ)'(B'B+D^2)^-1B)                   
    ∇_dlog(q_λ^a(θ)) =   diag(-(B'B+D^2)^-1D+(B'B+D^2)^-1(θ-μ)(θ-μ)'(B'B+D^2)^-1D) .

Fast calculation of these gradients can be undertaken using the Woodbury formula; see <cit.> for further details.

For VA1, λ^b=∅, so that
∇_λ^blog q_λ^b(u)=0. 
For VA2, if c_t=logω_t, then 

    log q_λ^b(u)=∑_t=1^T(1/2z_t^2 -c_t
     -(z_t-η_t)^2/2exp(2c_t)-log(b_t-a_t) ) ,

with derivatives ∇_c_tlog q_λ^b(u)=
exp(-2c_t)(z_t-η_t)^2-1 and 
∇_η_tlog q_λ^b(u)=(z_t-η_t)/ω^2_t. Last,
for VA3, if Ω^-1=LL', then

    log q_λ^b(u)=log |L| -1/2(z-η)'LL'(z-η)+∑_t=1^T 1/2z_t^2-log(b_t-a_t) ,

with derivatives 
∇_ηlog q_λ^b(u)=(z-η)'LL', and 

    ∇_(L)log q_λ^b(u)=
    ((L^-1)')'-1/2( (z-η)'⊗ (z-η))(I_T^2+K_T,T)
    (L⊗ I_T)

where K_T,T is a commutation matrix. Note that the gradient is for a full factor L, although
for the sparse L employed here we compute ∇_(L)log q_λ^b(u)
 using sparse matrix operations in Matlab, and only evaluate it
for the non-zero elements of L. Last, to derive this derivative we have used the identity ∂/∂(A)|A|=(|A|(A^-1)')' for invertible square matrix A.


14pt 

apa
