*****************************************
** Dynamic Conditional Correlations    **
**     in Political Science            **
** Matt Lebo and Jan Box-Steffensmeier **
** AJPS 52(3)  July, 2008              **
*****************************************
* Email matthew.lebo@sunysb.edu for questions
* NOTE: WITH EACH NEW VERSION OF RATS THERE ARE CHANGES TO THE
* UNDERLYING ESTIMATION OF GARCH MODELS, ESPECIALLY THE METHOD OF LIKELIHOOD
* ESTIMATION AND THE CALCULATION OF STANDARD ERRORS.
* YOU MIGHT GET RESULTS SLIGHTLY DIFFERENT FROM THE ARTICLE DUE
* TO THESE CHANGES.
* This file will run ALL analyses and figures for the paper, but the data are in separate files.

*********************************************************************************
* This file begins by making the first 2 figures, the data analyses follow that.*
*********************************************************************************
* TO CREATE FIGURES 1 and 2
* SIMULATE TWO SERIES WITH NO RELATIONSHIP
end 1
all 250
** create a random error component -- here distributed normally with mean = 0 and variance = 1
* This is the seed value (starting value) for the random number generator for the exact series in the paper
seed 11118
set u = %ran(1)
set v = %ran(.5)
** Set an intial value for y
set y1 = 0
set y2 = 0
** set a value for y in subsequent periods
set y1 2 * = u
set y2 2 * = v
* Create the new versions of the series, with volatile periods
SET y1vol = y1
SET y1vol 100 150 = 5*y1
SET y2vol = y2
SET y2vol 100 150 = 5*y2
* To make Figure 1:
set voltime = t>=100.and.t<=150
SPGRAPH(vfields=2,hfields=1,HEADER='Figure 1: Simulated Series Without and With Volatility')
graph(header='Two Simulated Series with Correlation = 0',VLABEL='Value of Series',max=10,min=-12.5,HLABEL='Time Point',klabels=||'Series 1               .  ','Series 2'||,key=loleft,patterns) 2
# y1
# y2
graph(shading=voltime,header='The Same Series with a (shaded) Period of High Volatility', $
 VLABEL='Value of Series',HLABEL='Time Point',klabels=||'Series 1 + Volatility','Series 2 + Volatility'||,key=loleft,patterns) 2
# y1vol
# y2vol
spgraph(done)

* Moving window estimation on series without volatility
clear coef1
clear tstat
do regend=51,250
 linreg y1 regend-50 regend
 # constant y2
 compute coef1(regend) = %beta(2)
 compute tstat(regend) = %tstats(2)
end do
print / coef1
* Moving window estimation on series with volatility
clear coef1vol
clear tstatvol
do regend=51,250
 linreg y1vol regend-50 regend
 # constant y2vol
 compute coef1vol(regend) = %beta(2)
 compute tstatvol(regend) = %tstats(2)
end do
print / tstatvol

* Kalman Filter estimation on series without volatility
equation kalcoef1 y1
# constant y2
system kalcoef1
end(system)
estimate(cohistory=cokalman1) 1 25
do entry=26, 250
   kalman(cohistory=cokalman1)
end do entry

* Kalman Filter estimation on series with volatility
equation kalcoef2 y1vol
# constant y2vol
system kalcoef2
end(system)
estimate(cohistory=cokalman2) 1 25
do entry=26, 250
   kalman(cohistory=cokalman2)  3
end do entry

* This makes Figure 2.  To make any single panel, just use those particular lines beginning with "graph"
SPGRAPH(vfields=3,hfields=1,HEADER='Figure 2: Moving Window and Kalman Filter Estimates')
graph(shading=voltime,header='Moving Window Estimates of the Regression Coefficient',subheader='Original Series and Volatile Series (over shaded period)', $
 VLABEL='Regression Coefficient',klabels=||'MW-50: Original','MW-50: Volatile'||,key=upright,patterns) 2
# coef1 50 250
# coef1vol 50 250
graph(shading=voltime,header='Moving Window Regression T-Statistics', $
 VLABEL='T-Statistic',klabels=||'MW-50: Original','MW-50: Volatile'||,key=upleft,patterns) 2
# tstat 50 250
# tstatvol 50 250
graph(shading=voltime,header='Kalman Filter Estimates of the Regression Coefficient', $
 VLABEL='Kalman Coefficient',HLABEL='Time Period',klabels=||'KF: Original','KF: Volatile'||,key=upright,patterns) 2
# cokalman1(2) 50 250
# cokalman2(2) 50 250
spgraph(done)
*********************************************************
* Next, the economic voting example         			     *
*********************************************************
* Beginning with the Tse Test - Table 1
********************************************************************************************
* Tse LM test for constant correlation.
* Tse, Y. K. (2000), "A Test for Constant Correlations in a Multivariate GARCH Model",
* Journal of Econometrics 98, 107-127.
* This is based on ESTIMA's source file.
* Revision Schedule:
*   07/05 Rewritten to use GARCH instruction
*
end 1
calendar 1978 1 12
allocate 0 2004:07
open data pres7804.xls
data(format=xls,org=obs)
* The following source files are needed to get estimates of the fractinal differencing parameter
* They are available at Estima.com
source(noecho) rgser.src
source(noecho) fif.src
* Get the fractionally differenced, stationary, versions of the variables before testing
difference ICS / ICSd
@rgser ICSd
@fif(d=-.11) ICSd / ICSdf
difference Npros5 / Npros5d
@rgser Npros5d
@fif(d=-.25) Npros5d / Npros5df
difference Nret / Nretd
@rgser Nretd
@fif(d=.10) Nretd / Nretdf
difference Ppros / Pprosd
@rgser Pprosd
@fif(d=-.33) Pprosd / Pprosdf
difference Presap / Presapd
@rgser Presapd
@fif(d=-.05) Presapd / Presapdf
difference Pret / Pretd
@rgser Pretd
@fif(d=-.24) Pretd / Pretdf
* Here's the test
compute n=2
*
* This has to follow estimation of a constant correlation model. The GARCH instruction
* below returns the derivatives of the log likelihood with respect to the CC model parameters
* as the vector[series] ccderives. We now need the derivatives with respect to the "delta"
* parameters in the extended model rho(i,j)(t)=rho(i,j)+delta(i,j)*e(i)(t)e(j)(t) at delta=0,
* which will be the CC estimates.
*
* For the Tse test of the correlation between presapdf and icsdf (Table 1, line 1) use this line:
garch(p=2,q=1,asymmetric,nomean,mv=cc,derives=ccderives,rvector=uv,hmat=h) / presapdf icsdf
* For the Tse test of the correlation between presapdf and 5-year national prospections (Table 1, line 2) use this line:
*garch(p=2,q=1,asymmetric,nomean,mv=cc,derives=ccderives,rvector=uv,hmat=h) / presapdf npros5df
* For the Tse test of the correlation between presapdf and national retrospections (Table 1, line 3) use this line:
*garch(p=2,q=1,asymmetric,nomean,mv=cc,derives=ccderives,rvector=uv,hmat=h) / presapdf nretdf
* For the Tse test of the correlation between presapdf and personal prospections (Table 1, line 3) use this line:
*garch(p=2,q=1,asymmetric,nomean,mv=cc,derives=ccderives,rvector=uv,hmat=h) / presapdf pprosdf
* For the Tse test of the correlation between presapdf and personal retrospections (Table 1, line 3) use this line:
*garch(p=2,q=1,asymmetric,nomean,mv=cc,derives=ccderives,rvector=uv,hmat=h) / presapdf pretdf
*
* This drops the first data point since the lagged ee' for that is just the initial
* estimate.
*
compute gstart=%regstart()+1,gend=%regend()
*
* This pulls the subdiagonal of the CC matrix out of the beta vector.
* The "0" in the nslot is equal to the number of mean model parameters
* that are estimated. This will be n if you estimate constant means.
*
dec symm qcinv(n,n)
compute nslot=0+3*n
ewise qcinv(i,j)=%if(i==j,1,%beta(nslot=nslot+1))
compute qcinv=inv(qcinv)
dec symm[series] xd(n-1,n-1)
*
dec vector ihu
do time=gstart+1,gend
   compute ux=uv(time),uux=%outerxx(uv(time-1)),hx=h(time),ihu=%diag(%sqrt(%xdiag(hx)))*inv(hx)*ux
   do i=1,n-1
      do j=1,i
         set xd(i,j) time time = -qcinv(i+1,j)*uux(i+1,j)+uux(i+1,j)*ihu(i+1)*ihu(j)
      end do j
   end do i
end do time
dec vector ihu
*
* The t-stat on the xd is the LM statistic -- divide by 2 for the p-value
*
set ones gstart gend = 1.0
linreg ones gstart+1 gend
# ccderives xd

***********************************************************
* Now the estimates of the DCC Models
* Just to be sure RATS' memory isn't muddled start from scratch.
end 1
calendar 1978 1 12
allocate 0 2004:07
open data pres7804.xls
data(format=xls,org=obs)
source(noecho) rgser.src
source(noecho) fif.src
* Make the series stationary prior to running the DCC models
difference ICS / ICSd
@rgser ICSd
@fif(d=-.11) ICSd / ICSdf
difference Npros5 / Npros5d
@rgser Npros5d
@fif(d=-.25) Npros5d / Npros5df
difference Nret / Nretd
@rgser Nretd
@fif(d=.10) Nretd / Nretdf
difference Ppros / Pprosd
@rgser Pprosd
@fif(d=-.33) Pprosd / Pprosdf
difference Presap / Presapd
@rgser Presapd
@fif(d=-.05) Presapd / Presapdf
difference Pret / Pretd
@rgser Pretd
@fif(d=-.24) Pretd / Pretdf

*GARCH MODELS
garch(p=2,q=1,asymetric,resides=r,iters=1000,hseries=h) / presapdf

********************************************************************
*** DCC MODEL FOR ICS and PRES APPROVAL -- FRACTIONALLY DIFFERENCED
********************************************************************
compute n=2
dec vect[series] x(n)
set x(1) = icsdf
set x(2) = Presapdf
*
dec vect[series] eps(n)
dec vect fullbeta(4*n+2)
* Do univariate GARCH models. Save the standardized residuals into eps(i).
* Copy the coefficients into the proper slots in the full beta matrix
*
do i=1,n
   garch(p=2,q=1,asymmetric,resides=r,hseries=h) / x(i)
   set eps(i) = r/sqrt(h)
   do j=1,4
      compute fullbeta(n*(j-1)+i)=%beta(j)
   end do j
end do i
* Note that "D" in the output is the asymetric parameter (m) and it must be multiplied by -1.
* compute the covariance matrix of the standardized residuals
* This gives the overall correlation between the series, R-bar, given in Table 1 and bottom of Table 2a.  Bottom left of matrix.
vcv(matrix=rr)
# eps
*
* Create the series[symm] uu (outer product of the residuals).
* Make it the unconditional value prior to the sample.
*
dec series[symm] uu q
gset uu %regstart() %regend() = %outerxx(%xt(eps,t))
gset uu 1 %regstart()-1 = rr
gset q = rr
*
* Log likelihood for the DCC phase, taking the residuals as given
* This gives DCC alpha and beta estimates
nonlin a b
dec frml[symm] qf
frml qf   = (qx=(1-a-b)*rr+a*uu{1}+b*q{1})
frml logl = q=qf,%logdensity(%cvtocorr(q),%xt(eps,t))
compute b=.80,a=.10
maximize logl 2 *
*
* Compute the estimates into the final two slots in fullbeta
*
compute fullbeta(4*n+1)=%beta(1),fullbeta(4*n+2)=%beta(2)
*
* Do one iteration of the full model with METHOD+BHHH to get the grand covariance matrix
* This is the one-step estimation
garch(p=1,q=1,mv=dcc,method=simplex,initial=fullbeta,iters=1000) / x
* Pull out the correlations
set dcorr 2 * = q{0}(1,2)/sqrt(q{0}(1,1)*q{0}(2,2))
* Figure 1 with just the top panel
set dempres = t>=1978:01.and.t<=1981:01.or.t>=1989:01.and.t<=1992:12.or.t>=2001:01
graph(shading=dempres,VLABEL='Correlation',HLABEL='Year',header='Figure 1: Dynamic Correlations', $
 subhead='Presidential Approval and Index of Consumer Sentiments') 1
# dcorr

* Kalman estimates of the same relationship
*
equation kalcoef1 presapdf
# constant icsdf
system kalcoef1
end(system)
estimate(cohistory=cokalman2) 1 5
do entry=6, 319
   kalman(cohistory=cokalman2)
end do entry
* Figure 3 with both panels
SPGRAPH(vfields=2,hfields=1,HEADER='Figure 3: Presidential Approval and the ICS, 1978-2004')
graph(shading=dempres,header='Dynamic Correlations',subheader='Shading separates administrations',vlabel='Dynamic Correlations') 1
# dcorr
graph(shading=dempres,header='Kalman Filter Coefficients',vlabel='Kalman Coefficient',HLABEL='Year') 1
# cokalman2(2)
spgraph(done)

** Next, make Table 3
***************************************************************************
*** DCC MODEL FOR National Prospections and PRES APPROVAL -- F. DIFFERENCED
***************************************************************************
compute n=2
dec vect[series] x(n)
set x(1) = npros5df
set x(2) = Presapdf
*
dec vect[series] eps(n)
dec vect fullbeta(4*n+2)
* Do univariate GARCH models. Save the standardized residuals into eps(i).
* Copy the coefficients into the proper slots in the full beta matrix
*
do i=1,n
   garch(p=2,q=1,asymmetric,resides=r,hseries=h) / x(i)
   set eps(i) = r/sqrt(h)
   do j=1,4
      compute fullbeta(n*(j-1)+i)=%beta(j)
   end do j
end do i
*
* compute the covariance matrix of the standardized residuals
*
vcv(matrix=rr)
# eps
*
* Create the series[symm] uu (outer product of the residuals).
* Make it the unconditional value prior to the sample.
*
dec series[symm] uu q
gset uu %regstart() %regend() = %outerxx(%xt(eps,t))
gset uu 1 %regstart()-1 = rr
gset q = rr
*
* Log likelihood for the DCC phase, taking the residuals as given
*
nonlin a b
dec frml[symm] qf
frml qf   = (qx=(1-a-b)*rr+a*uu{1}+b*q{1})
frml logl = q=qf,%logdensity(%cvtocorr(q),%xt(eps,t))
compute b=.80,a=.10
maximize logl 2 *
*
* Compute the estimates into the final two slots in fullbeta
*
compute fullbeta(4*n+1)=%beta(1),fullbeta(4*n+2)=%beta(2)
*
* Do one iteration of the full model with METHOD+BHHH to get the grand covariance matrix
*
garch(p=1,q=1,mv=dcc,method=bhhh,initial=fullbeta,iters=1000) / x

set dcorr 2 * = q{0}(1,2)/sqrt(q{0}(1,1)*q{0}(2,2))
* Make Figure 4
set dempres = t>=1978:01.and.t<=1981:01.or.t>=1989:01.and.t<=1992:12.or.t>=2001:01
SPGRAPH(vfields=2,hfields=1,HEADER='Figure 4: The Varying Efficacy of National Prospections, 1978-2004')
graph(shading=dempres,header='National Prospections and Presidential Approval',VLABEL='Approval Rating', $
 OVLABEL='Prospections Index',HLABEL='Year', $
 klabels=||'Presidential Approval','National Prospections'||,key=loright,patterns,style=line,overlay=line) 2
# presap
# npros5
graph(shading=dempres,header='Dynamic Correlations',VLABEL='Correlation',HLABEL='Year') 1
# dcorr
spgraph(done)

*****************************************************
** NOW THE IR EXAMPLE 										**
*****************************************************
* Clear and read in the new data.
end 1
calendar 1979 4 12
allocate 0 2004:06
open data levant.79-04.06.xls
data(format=xls,org=obs)
* Date ISRPAL ISRJOR ISRLEB ISRSYR ISRUAR ISRUSA ISRUSR *
*      PALISR PALJOR PALLEB PALSYR PALUAR PALUSA PALUSR *
*      JORISR JORPAL JORLEB JORSYR JORUAR JORUSA JORUSR *
*      LEBISR LEBPAL LEBJOR LEBSYR LEBUAR LEBUSA LEBUSR *
*      SYRISR SYRPAL SYRJOR SYRLEB SYRUAR SYRUSA SYRUSR *
*      UARISR UARPAL UARJOR UARLEB UARSYR UARUSA UARUSR *
*      USAISR USAPAL USAJOR USALEB USASYR USAUAR USAUSR *
*      USRISR USRPAL USRJOR USRLEB USRSYR USRUAR USRUSA *

* Read in the source files
source(noecho) rgser.src
source(noecho) fif.src

difference JORPAL / JORPALd
@rgser JORPALd
@fif(d=-.79) JORPALd / JORPALdf
difference PALJOR / PALJORd
@rgser PALJORd
@fif(d=-.61) PALJORd / PALJORdf
********************************************************************
*** DCC MODEL FOR PALJORDF and JORPALDF -- FRACTIONALLY DIFFERENCED
********************************************************************
compute n=2
dec vect[series] x(n)
set x(1) = paljordf
set x(2) = jorpaldf
*
dec vect[series] eps(n)
dec vect fullbeta(4*n+2)
* Do univariate GARCH models. Save the standardized residuals into eps(i).
* Copy the coefficients into the proper slots in the full beta matrix
*
do i=1,n
   garch(p=1,q=1,resides=r,method=bfgh,iterations=5000,piters=1000,hseries=h) / x(i)
   set eps(i) = r/sqrt(h)
   do j=1,4
      compute fullbeta(n*(j-1)+i)=%beta(j)
   end do j
end do i
*
* compute the covariance matrix of the standardized residuals
*
vcv(matrix=rr)
# eps
*
* Create the series[symm] uu (outer product of the residuals).
* Make it the unconditional value prior to the sample.
*
dec series[symm] uu q
gset uu %regstart() %regend() = %outerxx(%xt(eps,t))
gset uu 1 %regstart()-1 = rr
gset q = rr
*
* Log likelihood for the DCC phase, taking the residuals as given
*
nonlin a b
dec frml[symm] qf
frml qf   = (qx=(1-a-b)*rr+a*uu{1}+b*q{1})
frml logl = q=qf,%logdensity(%cvtocorr(q),%xt(eps,t))
compute b=.80,a=.10
maximize logl 2 *
*
* Compute the estimates into the final two slots in fullbeta
*
compute fullbeta(4*n+1)=%beta(1),fullbeta(4*n+2)=%beta(2)
*
* Do one iteration of the full model with METHOD+BHHH to get the grand covariance matrix
* If BHHH doesn't converge, SIMPLEX will
garch(p=1,q=1,mv=dcc,method=simplex,initial=fullbeta,iters=1000) / x
* Get the dynamic correlations
set dcorr 2 * = q{0}(1,2)/sqrt(q{0}(1,1)*q{0}(2,2))
print / dcorr
* Make Figure 5
graph(VLABEL='Correlation',HLABEL='Year',header='Figure 5: Dynamic Correlations', $
 subhead='Palestinian and Jordanian Interaction') 1
# dcorr