WIXLIB

forecasts of GDP growth have lower quantiles that decrease during recessionary periods, whereasthe upper quantiles remain stable over time, which is consistent with the findings of Adrian et al.(2017). The opposite is found for inflation, where the upper quantiles vary more over time than thelower quantiles, representing asymmetry in inflation uncertainty. Such features cannot be capturedby VARSV and UCSV models with symmetric disturbances. The real-time forecasts for the foureconomic variables from our copula model are compared to those from a VAR, VARSV, UCSV anda Gaussian copula model as in Smith and Vahey (2016). The point and density forecasts from ourcopula model are more accurate for inflation and the interest rate (the two most skewed variables)than those from the benchmark models. Notably, the current quarter forecasts of GDP growth are themost accurate, including when compared to those from the VARSV model, which is a high-performingforecasting model for this variable (Clark and Ravazzolo, 2015). We also construct bivariate realtime forecast densities for the change in the unemployment rate and GDP growth, which exhibitasymmetries and non-linearities. In line with Okun’s rule, Spearman correlations between the twovariables are negative, but show a severe weakening of the relationship during the Great Recession.Our model is not the only time series model that can produce asymmetric density forecasts. Forexample, a zero lower bound on the interest rate can be captured by adopting a shadow rate (Bauerand Rudebusch, 2016, Johannsen and Mertens, 2016), producing positively skewed density forecastsfor this variable. Another approach is to extend existing time series models using finite or infinitemixture models (Kalli and Griffin, 2018, Jensen and Maheu, 2010, Chiu et al., 2016). However, anadvantage of the copula model is that the margins can be modeled directly in an arbitrary mannerat an initial stage, and the dependence separately in a second stage.The rest of the paper is organized as follows. The heteroskedastic inversion copula is presented inSection 2, while Section 3 outlines the approach for its estimation and the computation of filtered andpredictive distributions. Section 4 presents the empirical analysis. This includes the data, existingevidence of asymmetry in these series, the results of fitting the copula model to the 2016Q2 vintage,and the real-time forecast study. Section 5 discusses implications and extensions.4

2Copula Model2.1Copula Time Series ModelConsider a continuous-valued multivariate time series {Yt }, with Yt (Y1,t , . . . , Yr,t )0 having r elements. Then Sklar’s theorem (Nelsen, 2006, p.46) states that the joint distribution function FY ofY (Y10 , . . . , YT0 )0 can be written asFY (y) C(u) ,(1)where y (y10 , . . . , yT0 )0 , yt (y1,t , . . . , yr,t )0 , u (u01 , . . . , u0T )0 , ut (u1,t , . . . , ur,t )0 , ui,t FYi,t (yi,t ),and FYi,t is the marginal distribution function of Yi,t . The function C is called the copula function,and is itself a joint distribution function on [0, 1]rT with strictly uniform margins (Nelsen, 2006, p.47).It captures the entire dependence structure of Y , which here is both the cross-sectional and temporaldependence structure of the time series. The major advantage of this decomposition is that it allowsthe margins FYi,t and dependence structure to be modeled separately.Throughout this paper we assume (strong) stationarity, so that the margins of the series aretime-invariant, and FYi,t Gi for i 1, . . . , r. In this case, the density of the data isfY (y) where gi (y) dG (y),dy ir YTY rTgi (yi,t ) ,FY (y) c(u) y1,1 · · · yr,Ti 1 t 1and c(u) rT u1,1 ··· ur,T(2)C(u) is called the ‘copula density’. The main challengeof this copula model is the selection of an appropriate copula function C (or, equivalently, densityc). For a univariate time series, Smith and Maneesoonthorn (2018) suggest constructing a copula byinverting a state space model, while for a multivariate time series Smith and Vahey (2016) construct acopula by inverting a VAR. We extend these ideas below, and propose using a new copula constructedfrom the inversion of a multivariate UCSV model.2.2Heteroskedastic Inversion CopulaConsider a second time series {Zt }, with Zt (Z1,t , . . . , Zr,t )0 having r elements. This follows astationary unobserved component model with stochastic volatility disturbances (UCSV) given by:Zt µt Λ0.5t εt ,5

µt pXBk µt k t ,k 1hi,t h̄i ai (hi,t 1 h̄i ) ηi,t , for i 1, . . . , r ,(3) Λt diag eh1,t , . . . , ehr,t , εt N (0, Ir ), ηt (η1,t , . . . , ηr,t )0 N (0, Ση ), and t N (0, Σ ). In thismodel, Zt is decomposed into the sum of two zero mean stationary components: a level componentµt which follows a VAR(p) model, and a variance component which has r correlated elements thateach follow stochastic volatility models, with Zi,t µi,t , hi,t N (µi,t , ehi,t ). Eq. (3) defines a nonlinearstate space model, with state vector (µ0t , h0t )0 , where ht (h1,t , . . . , hr,t )0 is the log-volatility vector.The autoregressive coefficients a (a1 , . . . , ar )0 are bounded so that 1 ai 1. We refer to {Zt }as a ‘pseudo time series’, because it is not observed directly. We extract its dependence structure byconstructing the T r-dimensional copula density of Z (Z10 , . . . , ZT0 )0 , and use it for c in Eq. (2).Let Ψ be the parameters of the UCSV model at Eq. (3), and Z have distribution function FZ (z Ψ).Because the pseudo time series is stationary, the marginal distribution function Fi (zi,t Ψ) of Zi,t istime invariant. Following Nelsen (2006, Sec.3.1), the copula function and density of this distributioncan be obtained by inversion of the usual expression of Sklar’s theorem, to give C(u Ψ) FZ F1 1 (u1,1 Ψ), F2 1 (u2,1 Ψ), . . . , Fr 1 (ur,T Ψ) Ψ , andfZ (z Ψ),QTi 1t 1 fi (zi,t Ψ)(4)c(u Ψ) Qrwhere zi,t Fi 1 (ui,t Ψ), zt (z1,t , . . . , zr,t )0 , z (z10 , . . . , zT0 )0 , fZ (z Ψ) fi (zi,t Ψ) dF (z Ψ).dzi,t i i,t rT z1,1 ··· zr,TFZ (z Ψ) andWe call this an ‘inversion copula’ (it is also called an ‘implicit copula’ byMcNeil et al. (2005, p.190)). However, for the data distribution of the UCSV model, fZ and fi cannotbe expressed in closed form, so that neither can c in Eq. (4). Nevertheless, the likelihood at Eq. (2)can still be computed by expressing it conditional on µ (µ01 , . . . , µ0T )0 and h (h01 , . . . , h0T )0 , andintegrating µ and h out using Bayesian methods, which is the approach we employ.The conditional likelihood can be obtained by a change of variables from Z to Y , as!TrYYgi (yi,t )f (y µ, h, Ψ) f (z µ, h, Ψ)JZ Y φr (zt ; µt , Λt ),(5)f (z Ψ)t 1i 1 i i,tQ Qgi (yi,t )where JZ Y Tt 1 ri 1 fi (zis the the Jacobian of the transformation, and φr (zt ; µt , Λt ) is thei,t Ψ)6

density of a N (µt , Λt ) distribution evaluated at zt . To compute Eq. (5), the most difficult task isevaluating the marginal densities fi , and the corresponding quantile functions Fi 1 . The latter isneeded to compute zi,t Fi 1 (ui,t Ψ) from ui,t Gi (yi,t ). The densityZfi (zi,t Ψ) φ1 (zi,t ; µi,t , hi,t )f (µi,t Ψ)f (hi,t Ψ)d(µi,t , hi,t ) .(6)The state µi,t Ψ N (0, di ), where d (d1 , . . . , dr ) can be computed from {B1 , . . . , Bp , Σ } (Lütkepohl,2005, pp.28-29), although (as we discuss below) we parameterize the vector autoregression directly22in terms of d instead. The state hi,t Ψ N (h̄i , s2i ), where s2i σ̄i,i/(1 a2i ), with σ̄i,ithe ith leadingdiagonal element of Ση . Substituting these two state densities into Eq. (6), and integrating out µi,tanalytically as a Gaussian, gives the expressionsZ fi (zi,t Ψ) φ1 zi,t ; 0, di eh φ1 h; h̄i , s2i dhZ Fi (zi,t Ψ) Φ1 zi,t ; 0, di eh φ1 h; h̄i , s2i dh ,(7)where Φr (z; m, Ω) is the distribution function of a N (m, Ω) evaluated at point z, and the integralsover h are computed numerically. Finally, to evaluate Fi 1 from Fi , we use the fast numerical methodproposed by Smith and Maneesoonthorn (2018, Appendix A).Tab. A1 in the Online Appendix depicts the transformations underlying the inversion copula timeseries model, and tabulates the various distribution and density functions.2.3Parameterization and IdentificationLet µ̃ (µ0p 1 , . . . , µ01 )0 , then the stationary VAR(p) model for the level component {µt } can beparameterized in terms of d and the unique partial correlations of the covariance matrix Var(µ̃). Ifthe partial correlation between series i and j, at lag k, is defined as0γi,j Corr (µi,t , µj,t intervening elements of vector µ̃) , if k 0 , andkγi,j Corr (µi,t k , µj,t intervening elements of vector µ̃) , for 1 k p ,kthen γ {γi,j; (i, j, k) P}, for P {1 j i r ; k 0} {1 i, j r; 1 k p}, isthe set of all unique partial correlations. There is a one-to-one relationship between {B1 , . . . , Bp , Σ }and {γ, d}. Lütkepohl (2005, pp.29-31) outlines how to compute {d, γ} from {B1 , . . . , Bp , Σ }, while7

Appendix A shows how to undertake the reverse. In total there are r2 p r(r 1)/2 partial correlations,so that there are 110 in our empirical analysis where r 5 and p 4. Adopting this parameterizationkis unconstrained on [ 1, 1] conditional onsimplifies our MCMC scheme in two ways. First, each γi,jthe other partials, simplifying generation and improving mixing. Second, the conditional posteriorof γ is not a function of F1 , . . . , Fr , so that generation of γ does not require re-computation ofz1 , . . . , zT , which is slow. Following Daniels and Pourahmadi (2009), we also parameterize the (r r)matrix Ση in terms of its variances and partial correlations.When forming a copula via inversion, all information about the marginal location and scale ofZi,t is lost. To identify the parameters, constraints on these moments are required; for example, inthe Gaussian copula Zi,t N (0, 1). We constrain the first two marginal moments of Zi,t . First,E(Zi,t ) 0 by definition of the model at Eq. (3). Second, the marginal variance isVar(Zi,t ) Var(µi,t ) Var(exp(hi,t /2) i,t ) di E(exp(hi,t )) di exp(h̄i s2i /2) ,2/(1 a2i ). Setting Var(Zi,t ) 1, and solving for h̄i , gives the equality constraintswhere s2i σ̄i,ih̄i log(1 di ) s2i /2, for i 1, . . . , r. These are simple to enforce in the MCMC scheme by evaluating(and not generating) h̄ (h̄1 , . . . , h̄r )0 at each sweep. In addition, because E(exp(hi,t )) 0, theparameters 0 di 1, for i 1, . . . , r. These inequality constaints are easily enforced in a Bayesianframework by adopting priors that constrain di (0, 1). With the parameterization and constraintsoutlined above, the copula parameters are now Ψ {γ, d, a, Ση }.It is important to note that the ‘role’ of each parameter in the UCSV model for {Zt } does nottranslate to the same role for {Yt }, so that Ψ is simply a set of copula parameters. However, di caninstead be interpreted as determining the relative importance of the two components for Zi,t : serialdependence is captured by the homoskedastic component µi,t when di 1, and the heteroskedasticcomponent when di 0. If all elements d 1, then E(exp(hi,t )) 0 and Zt µt . In this case, ourproposed copula reduces to the inversion copula of a Gaussian VAR(p), which is the Gaussian copuladiscussed by Smith and Vahey (2016), but without the D-vine representation employed by these1/2authors. If d 0, then µt 0 and Zt Λt εt , so that C is the inversion copula of a multivariate8

stochastic volatility model only.33.1Estimation and InferencePriorsAll elements of γ, d, a are assumed independent a priori, and constrained to their admissible domains.0Each cross-sectional partial correlation γi,jhas a N (0, 100) prior constrained to ( 1, 1), which is closekfor 1 k p has a N (0, 0.04to a flat prior. Each partial γi,j) prior constrained to ( 1, 1), so that thek2longer the lag of the partial correlation, the tighter is its prior density around zero. This is similar toa Minnesota prior (Litterman, 1986), where higher levels of shrinkage are applied to autoregressive2) I(0 s2i 3 ; ai coefficients of longer lags. The prior for di is uniform on [0, 1], while π(ai , σ̄i,i2} remain in a numerically stable region. Finally, uniform priors0.99) which ensures that {ai , σ̄i,ion ( 1, 1) are used for the r(r 1)/2 partial correlations of the matrix Ση , which is a standarduninformative prior on a correlation matrix (Daniels and Pourahmadi, 2009).3.2PosteriorEstimation of Ψ requires evaluation of the copula density at Eq. (4), which is the likelihood conditionalon the margins (G1 , . . . , Gr ). Even though the copula is high-dimensional, its parameters can stillbe estimated, as we now discuss. By exploiting the state space formulation of {Zt } in Eq (3), thenumerator of the copula density fZ (z Ψ) can be computed using a number of existing methods inthe state space literature. We adopt a Bayesian MCMC approach and explicitly generate the states{µ, h} as part of the following sampling scheme.Sampling SchemeStep 1: Generate from µ Ψ, h, y.Step 2: Generate from h Ψ, µ, y.Step 3: Generate from γ {Ψ\γ}µ, h, y.Step 4: Generate from di {Ψ\di }, µ, h, y, for i 1, . . . , r.Step 5: Generate from ai {Ψ\ai }, h, µ, y, for i 1, . . . , r.Step 6: Generate from Ση {Ψ\Ση }, h, µ, y.9

In Step 1, we compute {B1 , . . . , Bp , Σ } from {d, γ} as shown in the Appendix A, and then generateµ as a multivariate normal as in Chan and Jeliazkov (2009). To generate h in Step 2 we approximatethe distribution of log(ε2i,t ) for i 1, . . . , r, using a mixture of seven normals, and use the simulationsmoother. In Step 3, the partial correlations γ are sampled as pairs of elements chosen at random,where each pair is generated conditional on the remaining elements of γ using adaptive random walkMetropolis-Hastings (MH). In Steps 4, 5 and 6 the elements of d, a, and the diagonal elements of Σηare all generated one at a time using MH with normal approximations to the conditional posteriorsas the proposals. Finally, in Step 6 the partial correlations of the matrix Ση are generated oneat a time using adaptive random walk MH. After convergence, we collect the Monte Carlo sample{Ψ[1] , . . . , Ψ[J] } from which posterior inference is computed in the standard Bayesian fashion; theOnline Appendix gives a detailed description of the sampling scheme.3.3Simulation and ForecastingA major advantage of copulas constructed by inversion of a parametric model is that simulation of aniterate U C(· Ψ) is both fast and simple. To do so, first simulate a series of values Z1 , . . . , ZT fromthe UCSV model1 at Eq. (3), and then compute Ui,t Fi (Zi,t Ψ) for i 1, . . . , r and t 1, . . . , T .Monte Carlo estimates of a wide range of dependence measures — including slices of the copuladensity — can then be computed, as in Section 4.3.2.The h-step ahead posterior predictive density of YT h , conditional on Y1:T y, isZf (yT h y) f (yT h y, µ1:T , hT , Ψ)f (µ1:T , hT , Ψ y)d(µ1:T , hT , Ψ) .The integration above cannot be computed analytically, but it is fast and simple to simulate fromthe distribution as follows. At each sweep of the sampler, simulate from the conditional distribution ZT h µ1:T , hT , Ψ by first simulating µT h and hT h using the transition equations, andthen generating ZT h from the measurement equation, both given at Eq (3). Finally, computeYi,T h G 1i (Fi (Zi,T h Ψ)) for i 1, . . . , r.1To do so, we are careful to sample a pre-period to remove the effect of initial values of the VAR(p) process {µt }.10

The analogous (in-sample) h-step-ahead predictive density of Yt h , for t 1, . . . , T h, isZf (yt h y1:t ) f (yt h µ1:t , ht , Ψ)f (µ1:t , ht , Ψ y)d(µ1:t , ht , Ψ) ,which can be computed by simulating from Zt h µ1:t , ht , Ψ and then transforming, as above. Forboth densities, by simulating the iterates at each sweep of the sampler, Ψ is integrated out withrespect to its posterior in the usual Bayesian fashion.4Empirical Analysis4.1U.S. Data and AsymmetryWe consider quarterly data on U.S. output growth, inflation, the unemployment rate, a nominalshort-term interest rate and stock market volatility between 1954Q1 and 2016Q2. Output growth( ln(GDP)t ) and inflation (Inflt ) are the quarterly differences in the logarithms of real GDP andthe GDP price deflator, respectively. These are observed in real time, and are obtained from theReal-Time Data Set for Macroeconomists (Croushore and Stark, 2001). The data are provided inquarterly vintages, where vintage T 1 contains observations up to quarter T , and reflects the lagin the national accounts. Each vintage also provides revisions for the most recent quarters, andsometimes also the entire series. The interest rate (IRt ) is the quarterly average of the monthlysecondary market rate of the 3-Month T-Bill, and the unemployment rate (URt ) is the quarterlyaverage of the monthly seasonally-adjusted civilian unemployment rate; both obtained from theFRED database. Stock market volatility (VXOt ) is the quarterly average VXO volatility indexsourced from the FRED database, and backfilled prior to 1986 using calibrated daily SP500 realizedvolatility as in Bloom (2009). We set Yt ( ln(GDP)t , Inflt , URt , IRt , VXOt )0 , and Fig. 1(f–j) plotsthe five series for vintage 2016Q2.The four macroeconomic variables were considered by Clark and Ravazzolo (2015), Smith andVahey (2016) and Krüger and Nolte (2015). Some vintages are missing early observations of outputand the price deflator, and for these vintages we follow Clark and Ravazzolo (2015) and work onlywith the available information. For example, for the 1992Q1 vintage the 1959Q1 observations ofoutput and the price deflator are the earliest values used. Adrian et al. (2017) use the National11

Financial Conditions Index produced by the Federal Reserve Bank of Chicago as a measure of financialconditions. However, this only extends back to 1971, so that we instead use the VXO index, whichAdrian et al. (2017) also found to be related to GDP growth.There is strong evidence for asymmetry in the marginal distributions of these five variables.Shapiro-Wilk and Jarque-Bera tests (unreported) reject normality for all variables and vintages,except for some early vintages of the interest rate and unemployment. Using the nonparametrictest of symmetry of Bai and Ng (2005), symmetry is rejected at the 5% level for GDP growthvintages between 1977Q4 and 1991Q4, inflation from vintage 1992Q2 onwards, and stock marketvolatility until 1987Q4 and between 2006Q1 and 2008Q4. However, this test can have low power,and even strong deviations from symmetry may be undetected. Fig. 1(a–e) plots histograms of thevariables from vintage 2016Q2, and strong positive skew in inflation and market volatility is visible.Unemployment, market volatility and the interest rate are bounded on the left, with the latter havinga zero lower bound which it effectively reached between 2008Q4 and 2015Q4. Fig. 1 also shows theadaptive kernel density estimates (KDEs) of Shimazaki and Shinomoto (2010), which are stronglyasymmetric and heavy-tailed. They approximately capture lower bounds on unemployment and theinterest rate, capture the smooth tails and asymmetry in the variables, and have multiple modes forunemployment and the interest rate. Except where otherwise stated, KDEs from each vintage areused for the margins G1 , . . . , G5 in our copula time series models.Predictive distributions of growth can also exhibit asymmetry. To illustrate, we show that thepredictive distributions for GNP/GDP reported in the Survey of Professional Forecasters (SPF) areoften asymmetric. These are an ensemble of the survey respondents’ predictive distributions for the“average annual over average annual” growth rate for both the same and next year; see Smith andV

and Patton, 2017). However, as far as we are aware, ours is the rst copula constructed by inversion of a multivariate nonlinear time series model. The copula model is applied to quarterly U.S. real-time data on GDP growth, in ation, unem-ployment, a short term interest rate and stock market volatility. Adaptive kernel density estimates