new;

/* GAUSS code for bivariate VARFIMA(1,d,1) for the paper of Electoral Studies. 
   Here is for nonstationary economic expectations in Table 2 of Tsay (2012) 
   05232012
   Note: both AR and MA matrices are diagonal. 
   Author: Wen-Jen Tsay, Institute of Economics, Academia Sinica, Taiwan
   */
   

library maxlik, pgraph;
/*
#include pgraph.ext;
graphset;
*/
#include maxlik.ext;
maxset;

t_begin=hsec;


/* set ARMA(p,q) */
r = 2; /* dimensionality */
p = 1;
q = 1;
rmax = maxc(p|q);




_max_gradtol=1e-6;
_max_Algorithm=2;
_max_CovPar=1;
_max_MaxIters=250;
__output = 0;



load data[80,20] = dgw1.dat;

capp = data[.,3];
newapp = data[.,4];
econexp = data[.,5];
nytavg = data[.,6];
kg = data[.,10];
hb = data[.,14];
vetoes = data[.,16];
override = data[.,17];
intrasum = data[.,18];
mbills=data[.,19];


fitcapp = capp;
fiteconexp = econexp;


msize = rows(fitcapp);



trued1 = .23;
trued2 = 0.3;
truemu1 = 6;
truemu2 = 4;
truephi11 = .1805;
truephi12 = .1717;
truephi21 = .518;
truephi22 = .259;
rhohat = 0.355;
truesig11 = 2;
truesig22 = 1;
truesig12 = rhohat;
truebeta1 = 0.3;
truegama1 = 0.8;/* for nonstationary setting */
truetheta1 = 0.23;
truetheta2 = 0.2;



sigmaini = (truesig11|truesig12)
          ~(truesig12|truesig22)
          
;


sigini = chol(sigmaini); /*upper triangular */




  b00 = ( trued1/(1-2*trued1))|
          1|1|1|
          truephi11|
		  truephi22|
          truemu1|truemu2|
         (trued2/(1-2*trued2))|
         truebeta1|
		 truetheta1|
		 truetheta2
	
                 
;





count_npd=0;  /* counter for non-positive definite error */

est={};
estbstar={};
insamplefit={};



failed_npd=0;

output on;
output file = varfima1d1.out reset;

__title= "Conditonal Maximum Likelihood Estimator  ";
 { bstar,f,g,cov,ret } = maxprt((maxlik(fiteconexp~fitcapp,0,&lnlk, b00 )));
rett = ret;


   sigma1est = (bstar[2,1]|0)~(bstar[3,1]|bstar[4,1])
   ;
   sigmaest =  sigma1est'sigma1est;

   est = est|(
              (bstar[1,1]/( 1+2*abs(bstar[1,1]) ))|
              (sigmaest[1,1])|
              (sigmaest[2,2])|
              (sigmaest[1,2])|
              (bstar[5,1])|(bstar[6,1])|
              (bstar[7,1])|(bstar[8,1])|
			  (bstar[9,1])|
              (bstar[10,1]/( 1+2*abs(bstar[10,1]) ))|
			  (bstar[11,1]/( 1+abs(bstar[11,1]) ))|
			  (bstar[12,1]/( 1+abs(bstar[12,1]) ))|
    		 	 rett
            )'; 



?;
output file = varfima1d1.out on;
output on;
screen off;
format /rd 8,4;

"bstar is " bstar';
?;
"parameter estimate for double checking the following table of testing results";
est;
?;
?;
output off;


t_used=hsec-t_begin;





output file = varfima1d1.out on;
output on;
screen on;

/* to compute the information matrix by transforming the estimate from the first stage estimation*/


b1star = bstar[1,1]/(1+2*abs(bstar[1,1]));
b10star = bstar[10,1]/(1+2*abs(bstar[10,1]));


b2star = bstar[2,1];
b3star = bstar[3,1];
b4star = bstar[4,1];


sigma11star = (b2star|0)~(b3star|b4star);
/* smarter, to have p.d., upper tri */
sigma1star =  sigma11star'sigma11star;
sig11star = sigma1star[1,1];
sig22star = sigma1star[2,2];
sig12star = sigma1star[1,2];
sig21star = sigma1star[2,1];
b5star = bstar[5,1];
b6star = bstar[6,1];
b7star = bstar[7,1];
b8star = bstar[8,1];
b9star = bstar[9,1];
b11star =  bstar[11,1]/(1+abs(bstar[11,1]));
b12star =  bstar[12,1]/(1+abs(bstar[12,1]));


bnewstar = 
 b1star|sig11star|sig22star|sig12star|b5star|
            b6star|b7star|b8star|b9star|b10star|
           b11star|b12star

		   ;


_max_MaxIters=0;
_max_CovPar=1;
{ bstar,f,g,cov,ret } = maxprt(maxlik(fiteconexp[.,1]~fitcapp[.,1],0,&lnlk2nd,bnewstar ));



format /rd 8,4;
?;
?;
"Time used                ";;t_used/100;
?;
"total sample size  " rows(fiteconexp);






output off;


end;

/* MLE */
proc lnlk(b,zz);
local data,d1,d2,d3,
a1a,a2a,l1a,
b1a,b2a,l1b,
c1a,c2a,l1c,
l1z,l1ztil,
z1a,z2a,l1z2a,
z1atil,z2atil,l1z2atil,
/*
d2d12,l1d12,
d2d13,l1d13,
d2d23,l1d23,
*/
v00test,v00tilt,
b4,b5,b6,b7,b8,b9,
b10,b11,b12,b13,b14,b15,b16,b17,b18,
mu1hat,mu2hat,mu3hat,
beta1hat,gama1hat,
x1data,x2data,x3data,
x1data1,x2data1,x3data1,
sigma11,sigma1,
theta1hat,theta2hat,
sig11,sig22,sig33,sig12,sig13,sig21,sig23,sig31,sig32,
size,lag,
covar11,covar12,covar21,covar22,covar13,covar23,covar31,covar32,covar33,
i,crov,
v00,v00til,del00,del00til,
phi0,phi0til,v0,v0til,del0,del0til,vv,kk,vvtil,dd,ddtil,j,
vvfinal,v0final,
ff,dev,out1,out,var001,var00,ll;


d1 = b[1,1]/(1+2*abs(b[1,1]));
d2 = b[10,1]/(1+2*abs(b[10,1]));


/* for covariance matrix */
b4 = b[2,1];
b5 = b[3,1];
b6 = b[4,1];


b10 = b[5,1];
b14 = b[6,1];



mu1hat = b[7,1];
mu2hat = b[8,1];

beta1hat = b[9,1];
/*
gama1hat = b[13,1];
*/

theta1hat =  b[11,1]/(1+abs(b[11,1]));
theta2hat =  b[12,1]/(1+abs(b[12,1]));

gama1hat = 1;








sigma11 = (b4|0)~(b5|b6);

/* smarter, to have p.d., upper tri */
sigma1 =  sigma11'sigma11;
sig11 = sigma1[1,1];
sig22 = sigma1[2,2];
sig12 = sigma1[1,2];
sig21 = sigma1[2,1];

/* column 1 is econ expectation, column 2 is congressional approval */

x1data = zz[.,2] - mu1hat - beta1hat*zz[.,1];
x1data = x1data[2:rows(x1data),1];
x2data = zz[2:rows(zz),1] - mu2hat - gama1hat*zz[1:rows(zz)-1,1]; 



x1data1 = zeros(rows(x1data),1);
x2data1 = zeros(rows(x2data),1);


x1data1[1,.] = x1data[1,1];
x2data1[1,.] = x2data[1,1];

i = 2; do until i>rows(x1data);
x1data1[i,1] = x1data[i,1] -b10*x1data[i-1,1] /*- b11*x2data[i-1,1]*/;
x2data1[i,1] = x2data[i,1] /*-b13*x1data[i-1,1]*/ - b14*x2data[i-1,1];

i=i+1;endo;




/* data = zz;*/
data = x1data1~x2data1;
size = rows(data);
lag = size - 1;


{covar11} = covar01(size,d1,d1,theta1hat,theta1hat);
covar11 = sig11*covar11;


{covar12} = covar01(size,d1,d2,theta1hat,theta2hat);
covar12 = sig12*covar12;



{covar21} = covar01(size,d2,d1,theta2hat,theta1hat);
covar21 = sig21*covar21;


{covar22} = covar01(size,d2,d2,theta2hat,theta2hat);
covar22 = sig22*covar22;





/* note that autocovariance number is from 0 to nobs-1 */
/* if r is changed, this part needs to be adjusted */

crov = zeros(r,r*size);
i = 1;do until i>size;
crov[.,1+r*(i-1):r*i] =
                        (covar11[i,1]~covar12[i,1])|
                        (covar21[i,1]~covar22[i,1]);
                        
i = i+1;endo;




v00 = crov[.,1:r];
v00til = v00;
del00 = crov[.,1+r:2*r];
del00til = del00;


phi0 = zeros(r*lag,r*lag);
phi0til = zeros(r*lag,r*lag);


/* from 0 to size - 2 */
v0 = zeros(r,r*lag);
v0til = zeros(r,r*lag);

del0 = zeros(r,r*lag);
del0til = zeros(r,r*lag);


/* 11.4.13, p.422 */
v0[.,1:r] = v00;
v0til[.,1:r] = v00til;
del0[.,1:r] = del00;
del0til[.,1:r] = del00til;



phi0[1:r,1:r] = del00*inv(v00til);
phi0til[1:r,1:r] = del00til*inv(v00);


i = 2;do until i>lag;

clear vv;
kk=1;do until kk>i-1;
vv = vv +
        phi0[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
        (crov[.,1+r*(kk):r*(kk+1)]');

kk = kk+1;endo;

v0[.,1+r*(i-1):r*i] = crov[.,1:r] - vv;


clear vvtil;
kk=1;do until kk>i-1;
vvtil = vvtil +
        phi0til[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
        (crov[.,1+r*(kk):r*(kk+1)]);

kk = kk+1;endo;

v0til[.,1+r*(i-1):r*i] = crov[.,1:r] - vvtil;




clear dd;
kk=1;do until kk>i-1;
dd = dd +
        phi0[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
       (crov[.,1+r*(i-kk):r*(i-kk+1)]);

kk = kk+1;endo;

del0[.,1+r*(i-1):r*i] = crov[.,1+r*(i):r*(i+1)] - dd;


clear ddtil;
kk=1;do until kk>i-1;
ddtil = ddtil +
        phi0til[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
       (crov[.,1+r*(i-kk):r*(i-kk+1)]');

kk = kk+1;endo;

del0til[.,1+r*(i-1):r*i] = crov[.,1+r*(i):r*(i+1)]' - ddtil;




/* 11.4.3 */

              phi0[1+r*(i-1):r*i,1+r*(i-1):r*i] =
                          del0[.,1+r*(i-1):r*i]
                    *inv(v0til[.,1+r*(i-1):r*i]);


           phi0til[1+r*(i-1):r*i,1+r*(i-1):r*i] =
                       del0til[.,1+r*(i-1):r*i]
                       *inv(v0[.,1+r*(i-1):r*i]);



j = 1;do until j>i-1;


              phi0[1+r*(j-1):r*j,1+r*(i-1):r*i] =
              phi0[1+r*(j-1):r*j,1+r*(i-1-1):r*(i-1)]-
              phi0[1+r*(i-1):r*i,1+r*(i-1):r*i]*
           phi0til[1+r*(i-j-1):r*(i-j),1+r*(i-1-1):r*(i-1)];

           phi0til[1+r*(j-1):r*j,1+r*(i-1):r*i] =
           phi0til[1+r*(j-1):r*j,1+r*(i-1-1):r*(i-1)]-
           phi0til[1+r*(i-1):r*i,1+r*(i-1):r*i]*
              phi0[1+r*(i-j-1):r*(i-j),1+r*(i-1-1):r*(i-1)];


j=j+1;endo;
i = i+1;endo;

/* the n-1: final v */
clear vvfinal;
kk=1;do until kk>lag;
vvfinal = vvfinal +
        phi0[1+r*(kk-1):r*kk,1+r*(lag-1):r*lag]*
        (crov[.,1+r*(kk):r*(kk+1)]');

kk = kk+1;endo;

v0final = crov[.,1:r] - vvfinal;




dev = zeros(r,size-1);
i = 1; do until i > size-1;

clear ff;
kk=1;do until kk>i;
ff = ff +
        phi0[1+r*(kk-1):r*kk,1+r*(i-1):r*i]*(data[i-kk+1,.]');

kk = kk+1;endo;

dev[.,i] = data[i+1,.]' - ff;

i = i+1;endo;

dev = (data[1,.]')~dev;


  out1 = zeros(size-1,1);
  i = 1;do until i > size-1;
  out1[i,1] = dev[.,i]' * inv(v0[.,1+r*(i-1):i*r]) * dev[.,i];
  i = i + 1;
  endo;

  out = out1|
  (dev[.,size]' * inv(v0final) * dev[.,size]);

  var001 = zeros(size-1,1);
  i = 1;do until i > size-1;
  var001[i,1] = det(v0[.,1+r*(i-1):i*r]);
  i = i + 1;
  endo;

  var00 = var001|(det(v0final));


  ll = -0.5*r*ln(2*pi) -0.5*ln( var00 ) - 0.5*out;



retp(ll);
endp;



fn covar(dm,dn,nobs) =
     (gamma(2-dm-dn)/( (1-dm-dn)*gamma(1-dm)* gamma(1-dn) ) )
/*     (gamma(1-dm-dn)/( gamma(1-dm)* gamma(1-dn) ) )*/

 .* (1|cumprodc(seqa(dn,1,nobs-1)./seqa(-dm+1,1,nobs-1)))
;




proc(1) = covar01(nobs,dm,dn,t1,t2);

local dmm,dnn, codevar, c1, d1, rj;

dmm = (dm + conj(dm))/2;
dnn = (dn + conj(dn))/2;


codevar =
     (gamma(2-dmm-dnn)/( (1-dmm-dnn)*gamma(1-dmm)* gamma(1-dnn) ) )
/*     (gamma(1-dm-dn)/( gamma(1-dm)* gamma(1-dn) ) )*/

 .* (1|cumprodc(seqa(dnn,1,nobs-1)./seqa(-dmm+1,1,nobs-1)))
;

c1 = seqa(-dmm,1,nobs)./seqa(dnn-1,1,nobs);
d1 = seqa(dnn,1,nobs)./seqa(1-dmm,1,nobs);


rj = (1+t1*t2) + (t1.*c1) + (t2 .* d1);
codevar = codevar.*rj;

retp(codevar);
endp;












/* MLE */
proc lnlk2nd(b,zz);
local data,d1,d2,d3,
a1a,a2a,l1a,
b1a,b2a,l1b,
c1a,c2a,l1c,
l1z,l1ztil,
z1a,z2a,l1z2a,
z1atil,z2atil,l1z2atil,
v00test,v00tilt,
b4,b5,b6,b7,b8,b9,
b10,b11,b12,b13,b14,b15,b16,b17,b18,
mu1hat,mu2hat,mu3hat,
beta1hat,gama1hat,
x1data,x2data,x3data,
x1data1,x2data1,x3data1,
sigma11,sigma1,
sig11,sig22,sig33,sig12,sig13,sig21,sig23,sig31,sig32,
size,lag,
theta1hat,theta2hat,
covar11,covar12,covar21,covar22,covar13,covar23,covar31,covar32,covar33,
i,crov,
v00,v00til,del00,del00til,
phi0,phi0til,v0,v0til,del0,del0til,vv,kk,vvtil,dd,ddtil,j,
vvfinal,v0final,
ff,dev,out1,out,var001,var00,ll;


d1 = b[1,1];
d2 = b[10,1];



b4 = b[2,1];
b5 = b[3,1];
b6 = b[4,1];


b10 = b[5,1];
b14 = b[6,1];



mu1hat = b[7,1];
mu2hat = b[8,1];

beta1hat = b[9,1];
/*
gama1hat = b[13,1];
*/
gama1hat = 1;

theta1hat = b[11,1];
theta2hat = b[12,1];






sig11 = b[2,1];
sig22 = b[3,1];
sig12 = b[4,1];
sig21 = sig12;

/* column 1 is econ expectation, column 2 is congressional approval */

x1data = zz[.,2] - mu1hat - beta1hat*zz[.,1];
x1data = x1data[2:rows(x1data),1];
x2data = zz[2:rows(zz),1] - mu2hat - gama1hat*zz[1:rows(zz)-1,1]; 



x1data1 = zeros(rows(x1data),1);
x2data1 = zeros(rows(x2data),1);


x1data1[1,.] = x1data[1,1];
x2data1[1,.] = x2data[1,1];

i = 2; do until i>rows(x1data);
x1data1[i,1] = x1data[i,1] -b10*x1data[i-1,1] /*- b11*x2data[i-1,1]*/;
x2data1[i,1] = x2data[i,1] /*-b13*x1data[i-1,1]*/ - b14*x2data[i-1,1];

i=i+1;endo;




/* data = zz;*/
data = x1data1~x2data1;
size = rows(data);
lag = size - 1;





{covar11} = covar01(size,d1,d1,theta1hat,theta1hat);
covar11 = sig11*covar11;


{covar12} = covar01(size,d1,d2,theta1hat,theta2hat);
covar12 = sig12*covar12;



{covar21} = covar01(size,d2,d1,theta2hat,theta1hat);
covar21 = sig21*covar21;



{covar22} = covar01(size,d2,d2,theta2hat,theta2hat);
covar22 = sig22*covar22;






/* note that autocovariance number is from 0 to nobs-1 */
/* here for different r, need to be adjusted */

crov = zeros(r,r*size);
i = 1;do until i>size;
crov[.,1+r*(i-1):r*i] =
                        (covar11[i,1]~covar12[i,1])|
                        (covar21[i,1]~covar22[i,1]);
                        
i = i+1;endo;







v00 = crov[.,1:r];
v00til = v00;
del00 = crov[.,1+r:2*r];
del00til = del00;


phi0 = zeros(r*lag,r*lag);
phi0til = zeros(r*lag,r*lag);


/* from 0 to size - 2 */
v0 = zeros(r,r*lag);
v0til = zeros(r,r*lag);

del0 = zeros(r,r*lag);
del0til = zeros(r,r*lag);


/* 11.4.13, p.422 */
v0[.,1:r] = v00;
v0til[.,1:r] = v00til;
del0[.,1:r] = del00;
del0til[.,1:r] = del00til;



phi0[1:r,1:r] = del00*inv(v00til);
phi0til[1:r,1:r] = del00til*inv(v00);


i = 2;do until i>lag;

clear vv;
kk=1;do until kk>i-1;
vv = vv +
        phi0[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
        (crov[.,1+r*(kk):r*(kk+1)]');

kk = kk+1;endo;

v0[.,1+r*(i-1):r*i] = crov[.,1:r] - vv;


clear vvtil;
kk=1;do until kk>i-1;
vvtil = vvtil +
        phi0til[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
        (crov[.,1+r*(kk):r*(kk+1)]);

kk = kk+1;endo;

v0til[.,1+r*(i-1):r*i] = crov[.,1:r] - vvtil;




clear dd;
kk=1;do until kk>i-1;
dd = dd +
        phi0[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
       (crov[.,1+r*(i-kk):r*(i-kk+1)]);

kk = kk+1;endo;

del0[.,1+r*(i-1):r*i] = crov[.,1+r*(i):r*(i+1)] - dd;


clear ddtil;
kk=1;do until kk>i-1;
ddtil = ddtil +
        phi0til[1+r*(kk-1):r*kk,1+r*(i-1-1):r*(i-1)]*
       (crov[.,1+r*(i-kk):r*(i-kk+1)]');

kk = kk+1;endo;

del0til[.,1+r*(i-1):r*i] = crov[.,1+r*(i):r*(i+1)]' - ddtil;




/* 11.4.3 */

              phi0[1+r*(i-1):r*i,1+r*(i-1):r*i] =
                          del0[.,1+r*(i-1):r*i]
                    *inv(v0til[.,1+r*(i-1):r*i]);


           phi0til[1+r*(i-1):r*i,1+r*(i-1):r*i] =
                       del0til[.,1+r*(i-1):r*i]
                       *inv(v0[.,1+r*(i-1):r*i]);



j = 1;do until j>i-1;


              phi0[1+r*(j-1):r*j,1+r*(i-1):r*i] =
              phi0[1+r*(j-1):r*j,1+r*(i-1-1):r*(i-1)]-
              phi0[1+r*(i-1):r*i,1+r*(i-1):r*i]*
           phi0til[1+r*(i-j-1):r*(i-j),1+r*(i-1-1):r*(i-1)];

           phi0til[1+r*(j-1):r*j,1+r*(i-1):r*i] =
           phi0til[1+r*(j-1):r*j,1+r*(i-1-1):r*(i-1)]-
           phi0til[1+r*(i-1):r*i,1+r*(i-1):r*i]*
              phi0[1+r*(i-j-1):r*(i-j),1+r*(i-1-1):r*(i-1)];


j=j+1;endo;
i = i+1;endo;

/* the n-1: final v */
clear vvfinal;
kk=1;do until kk>lag;
vvfinal = vvfinal +
        phi0[1+r*(kk-1):r*kk,1+r*(lag-1):r*lag]*
        (crov[.,1+r*(kk):r*(kk+1)]');

kk = kk+1;endo;

v0final = crov[.,1:r] - vvfinal;




dev = zeros(r,size-1);
i = 1; do until i > size-1;

clear ff;
kk=1;do until kk>i;
ff = ff +
        phi0[1+r*(kk-1):r*kk,1+r*(i-1):r*i]*(data[i-kk+1,.]');

kk = kk+1;endo;

dev[.,i] = data[i+1,.]' - ff;

i = i+1;endo;

dev = (data[1,.]')~dev;


  out1 = zeros(size-1,1);
  i = 1;do until i > size-1;
  out1[i,1] = dev[.,i]' * inv(v0[.,1+r*(i-1):i*r]) * dev[.,i];
  i = i + 1;
  endo;

  out = out1|
  (dev[.,size]' * inv(v0final) * dev[.,size]);

  var001 = zeros(size-1,1);
  i = 1;do until i > size-1;
  var001[i,1] = det(v0[.,1+r*(i-1):i*r]);
  i = i + 1;
  endo;

  var00 = var001|(det(v0final));


  ll = -0.5*r*ln(2*pi) -0.5*ln( var00 ) - 0.5*out;



retp(ll);
endp;