gausset;
library cmlmt,pgraph;
#include cmlmt.sdf

/* K1    - Lag number in autocorrelation test                               */
/* K2    - Lag number in heteroscedasticity test                            */
let matK1=1 2 3 4 5 6 7 8 9 10 20 30;
let matK2=1 2 3 4 5 6 7 8 9 10 20 30;
let Tboots=20000;
/* UNIVARIATE TESTS                                                         */
/* Normality tests take                        1  columns                   */
/* Autocorrelation tests take         rows(matK1) columns                   */
/* Heteroscedasticity tests take      rows(matK2) columns                   */
Result=zeros(1,1+rows(matK1)+rows(matK2));
                                                    
/* The data is loaded                                                       */
load z[273,1]=EURUSD.txt;
z=z[1:252,1];

#include ds.sdf
struct DS z0;
z0=reshape(dsCreate,3,1);
z0[1].dataMatrix=z[3:rows(z)];
z0[2].dataMatrix=z[2:(rows(z)-1)];
z0[3].dataMatrix=z[1:(rows(z)-2)];

apuria=seqa(1989,1/12,rows(z0[1].dataMatrix)); 
graphset;
ylabel("Original Series");          
xlabel("Time");             
if cols(z0[1].dataMatrix)>1;
  _pltype = seqa(1,1,cols(z0[1].dataMatrix));
endif;
xtics(1989,2010,4,1);
_ptek="WL2AR2Series.tkf";
title("Time series of the data");     
xy(apuria,z0[1].dataMatrix);
ret = tkf2eps("WL2AR2Series.tkf", "WL2AR2Series.eps");

phi1=(0.019|1.26|-0.27);
phi2=(0.09);
s1=0.065;
s2=0.0149;
beta=zeros(2,1);


struct PV x0;
x0 = pvCreate;
x0 = pvPacki(x0,phi1,"phi1",1);
x0 = pvPacki(x0,phi2,"phi2",2);
x0 = pvPacki(x0,s1,"s1",3);
x0 = pvPacki(x0,s2,"s2",4);
x0 = pvPacki(x0,beta,"beta",5);

initialvalues = pvGetParVector(x0);

/* Setting variable values for ESTIMATION                                   */
struct cmlmtControl ctl; 
ctl = cmlmtControlCreate;
ctl.title = "Wong&Li model: 2 regimes, 2 lags in each estimated using the cml-routine";
ctl.Algorithm = 3;
/*  = 0,   Modifed BFGS
**  = 1,   BFGS (Broyden, Fletcher, Goldfarb, Shanno) (default)                      
**  = 2,   DFP (Davidon, Fletcher, Powell)                                  
**  = 3,   NEWTON (Newton-Raphson) 	                                      
**  = 4,   BHHH 	                                      */
ctl.LineSearch = 1;
/*  = 0, augmented trust region method (requires constraints)
**  = 1, STEPBT (quadratic and cubic curve ¯t) (default)
**  = 2, Brent's method
**  = 3, BHHHStep
**  = 4, half
**  = 5, Wolfe's condition */
ctl.PrintIters=1;
ctl.maxIters=1000;
ctl.Active={1,1,1,1,1,1,1,1}; 
ctl.Bounds = {-100 100, -2 2, -0.99 0.99, -100 100,
0 100, 0 100, -100 100, -100 100};
ctl.DirTol = 1e-4;
screen on;
struct cmlmtResults out1;
out1 = CMLmt(&ll,x0,z0,ctl);
x0=pvPutParVector(x0,(pvGetParvector(out1.Par)));

/* Log-Likelihood function for a NewMixture_of_WL2AR2_models                  */
proc ll(struct PV x0, struct DS z0, ind); 
local phi1,phi2, s1,s2, p, f, T, h11, h12,apu,beta,lags,lags2;
struct modelResults mm;
phi1 = pvUnpack(x0,1);
phi2 = pvUnpack(x0,2);
phi2 = phi2|phi1[2:3];
s1 = pvUnpack(x0,3); 
s2 = pvUnpack(x0,4);
beta = pvUnpack(x0,5);
T=rows(z0[1].dataMatrix);
lags=ones(rows(z0[1].dataMatrix),1)~z0[2].dataMatrix~z0[3].dataMatrix;
lags2=ones(rows(z0[1].dataMatrix),1)~z0[3].dataMatrix;
h11=zeros(T,1); h12=h11;
for t (1,T,1);
  h11[t]=(2*pi*det(s1))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi1)'*invpd(s1)*(z0[1].dataMatrix[t]-lags[t,.]*phi1)/2);
  h12[t]=(2*pi*det(s2))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi2)'*invpd(s2)*(z0[1].dataMatrix[t]-lags[t,.]*phi2)/2);
endfor;    
apu=exp(lags2*beta);
p=apu./(ones(rows(apu),1)+apu);
local d,y1,y2,y3,n,muu1,muu2;
n=1;
d=zeros(T,n); y1=zeros(T,1); y2=zeros(T,n); y3=zeros(T,1);
d=z0[1].dataMatrix;
muu1=lags*phi1; muu2=lags*phi2;
for t (1,T,1);
  y1[t,1]=ln(g(d[t,.]',muu1[t],s1,muu2[t],s2,p[t]));
endfor;
if ind[1];
  mm.Function = y1;
endif;
retp(mm);
endp;

/* This procedure computes the gradient for the log-likelihood function     */
/* This procedure calculates the gradient for multivar. quantile residuals  */
/* This procedure calculates the gradient for joint quantile residuals      */
/* Returns a (T+nT+T) x 1 -matrix                                           */

proc fct1(struct PV x0, struct DS z0); 
local phi1,phi2, s1,s2, p, f, T, h11, h12,apu,beta,lags,lags2;
struct modelResults mm;
phi1 = pvUnpack(x0,1);
phi2 = pvUnpack(x0,2);
phi2 = phi2|phi1[2:3];
s1 = pvUnpack(x0,3); 
s2 = pvUnpack(x0,4);
beta = pvUnpack(x0,5);
T=rows(z0[1].dataMatrix);
lags=ones(rows(z0[1].dataMatrix),1)~z0[2].dataMatrix~z0[3].dataMatrix;
lags2=ones(rows(z0[1].dataMatrix),1)~z0[3].dataMatrix;
h11=zeros(T,1); h12=h11;
for t (1,T,1);
  h11[t]=(2*pi*det(s1))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi1)'*invpd(s1)*(z0[1].dataMatrix[t]-lags[t,.]*phi1)/2);
  h12[t]=(2*pi*det(s2))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi2)'*invpd(s2)*(z0[1].dataMatrix[t]-lags[t,.]*phi2)/2);
endfor;    
apu=exp(lags2*beta);
p=apu./(ones(rows(apu),1)+apu);
local d,y1,y2,n,muu1,muu2,be,apu2;
n=1;
d=zeros(T,n); y1=zeros(T,1); y2=zeros(T,n); 
d=z0[1].dataMatrix;
muu1=lags*phi1; muu2=lags*phi2;
for t (1,T,1);
  y1[t,1]=ln(g(d[t,.]',muu1[t]',s1,muu2[t]',s2,p[t]));
  y2[t,.]=resid(&g,&cumulg,d[t,.]',muu1[t],s1,muu2[t],s2,p[t])';
endfor; 
retp(y1|vecr(y2')); 
endp; 

proc fct2(struct PV x0, struct DS z0); 
local phi1,phi2, s1,s2, p, f, T, h11, h12,apu,beta,lags,lags2;
struct modelResults mm;
phi1 = pvUnpack(x0,1);
phi2 = pvUnpack(x0,2);
phi2 = phi2|phi1[2:3];
s1 = pvUnpack(x0,3); 
s2 = pvUnpack(x0,4);
beta = pvUnpack(x0,5);
T=rows(z0[1].dataMatrix);
lags=ones(rows(z0[1].dataMatrix),1)~z0[2].dataMatrix~z0[3].dataMatrix;
lags2=ones(rows(z0[1].dataMatrix),1)~z0[3].dataMatrix;
h11=zeros(T,1); h12=h11;
for t (1,T,1);
  h11[t]=(2*pi*det(s1))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi1)'*invpd(s1)*(z0[1].dataMatrix[t]-lags[t,.]*phi1)/2);
  h12[t]=(2*pi*det(s2))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi2)'*invpd(s2)*(z0[1].dataMatrix[t]-lags[t,.]*phi2)/2);
endfor;    
apu=exp(lags2*beta);
p=apu./(ones(rows(apu),1)+apu);
local d,y1,y2,y3,n,muu1,muu2;
n=1;
d=zeros(T,n); y1=zeros(T,1); y2=zeros(T,n); y3=zeros(T,1);
d=z0[1].dataMatrix;
muu1=lags*phi1; muu2=lags*phi2;
for t (1,T,1);
  y1[t,1]=ln(g(d[t,.]',muu1[t],s1,muu2[t],s2,p[t]));
endfor;
retp(y1); 
endp; 

proc datasimu_WL2AR2(struct PV x, Tboots); 
local phi1,phi2,gamma1,gamma2,s1,s2,beta;
phi1 = pvUnpack(x0,1);
phi2 = pvUnpack(x0,2);
s1 = pvUnpack(x0,3); 
s2 = pvUnpack(x0,4);
beta = pvUnpack(x0,5);
local epsilon, boolean, data, ey1,ey2,bool,apu,p;
epsilon=rndn(Tboots+205,1);
data=zeros(Tboots+205,1);
bool=zeros(Tboots+205,1);
apu=exp((1~0)*beta);
p=apu./(1+apu);
boolean=RNDBernoulli(p);
if boolean==1;
    data[1:2]=phi1[1]*ones(2,1)+sqrt(s1).*epsilon[1:2];
else;
    data[1:2]=phi2[1]*ones(2,1)+sqrt(s2).*epsilon[1:2];
endif;
apu=exp((1~data[1])*beta);
p=apu./(1+apu);
boolean=RNDBernoulli(p);
bool[2,1]=boolean;
for k (3,Tboots+205,1); 
  if boolean==1;
    data[k]=phi1[1]+phi1[2:3]'*data[k-1:k-2]+sqrt(s1).*epsilon[k];
  else;
    data[k]=phi2[1]+phi1[2:3]'*data[k-1:k-2]+sqrt(s2).*epsilon[k];
  endif;
  apu=exp((1~data[k-1])*beta);
  p=apu./(1+apu);
  boolean=RNDBernoulli(p);
  bool[k,1]=boolean;
endfor;
retp(data[201:Tboots+200]|bool[201:Tboots+200]);
endp;

proc RNDBernoulli(p);
local bool;
bool=rndu(1,1);
if p<0;   retp (error(0));
elseif p>1;   retp (error(0));
else;
  if bool<=p;    retp(1); 
  else;      retp(2);
  endif;
endif;
endp;

/* CALCULATION OF QUANTILE RESIDUALS                                        */
T=rows(z0[1].dataMatrix);
lqr=fct1(x0,z0);
mr=lqr[(T+1):2*T,.];

/* NUMERICAL DERIVATIVES                                                    */
#include optim.sdf 
deriv = gradMT(&fct1,x0,z0);
s = deriv[1:T,.];
mg = deriv[(T+1):2*T,.];
x0hessian = hessMT(&fct2,x0,z0);

/* BOOTSRAPPED QUANTILE RESIDUALS                                           */
struct DS z0boots;
z0boots=reshape(dsCreate,3,1);
apuri=datasimu_WL2AR2(x0,Tboots+2);
z0boots[1].dataMatrix=apuri[3:(Tboots+2)];
z0boots[2].dataMatrix=apuri[2:(Tboots+1)];
z0boots[3].dataMatrix=apuri[1:(Tboots)];
lqrboots=fct1(x0,z0boots);
mrboots=lqrboots[(Tboots+1):2*Tboots,.];

/* THEIR NUMERICAL DERIVATIVES                                              */
dim=rows(initialvalues);
n=cols(apuri);
sboots=zeros(Tboots,dim); mgboots=zeros(Tboots,n*dim);
derivboots = gradMT(&fct1,x0,z0boots);
sboots=derivboots[1:Tboots,.];
mgboots=derivboots[(Tboots+1):(2*Tboots),.];

/****************************************************************************/
/* UNIVARIATE NORMALITY TESTS                                               */ 
/****************************************************************************/
   
Result[1,1]= ntestingboots(mr,mrboots,sboots,mgboots);
  
/****************************************************************************/
/* UNIVARIATE AUTOCORRELATION TEST                                          */ 
/****************************************************************************/

Result[1,2:1+rows(matK1)]
=atestingboots(mr,mrboots,sboots,mgboots,matK1);
      
/***************************************************************************/
/* UNIVARIATE HETEROSCEDASTICITY TEST                                       */ 
/****************************************************************************/

Result[1,2+rows(matK1):1+rows(matK1)+rows(matK2)]
=htestingboots(mr,mrboots,sboots,mgboots,matK2);

let pr=1; 
outwidth 80;
output file = estim_Wong_Li_2Regim_AR2.asc on;
screen off; 

call cmlmtPrt(out1);
print "";
print "Initial values used in estimation";
print pvGetParNames(x0)'; 
print initialvalues';
print "";
print "Sample size used for estimation " ftos(T,"%*.*lf",1,0);
print "";
print "The value of the likelihood at minimum " ftos(out1.Fct,"%*.*lf",1,3);
print "";
print "The AIC value " ftos(-2*(out1.Fct)+2*(rows(initialvalues)-1),"%*.*lf",1,3);
print "";
print "The BIC value " ftos(-2*(out1.Fct)+(rows(initialvalues)-1)*
                         ln(T),"%*.*lf",1,3);
                                         
xy0=pvGetParVector(x0);
h11=(2*pi*xy0[7])^-0.5
     *exp(-(z0[1].dataMatrix-xy0[1]-xy0[2]*z0[2].dataMatrix
     -xy0[3]*z0[3].dataMatrix)^2/xy0[7]/2);
h12=(2*pi*xy0[8])^-0.5
     *exp(-(z0[1].dataMatrix-xy0[4]-xy0[5]*z0[2].dataMatrix
     -xy0[6]*z0[3].dataMatrix)^2/xy0[8]/2);

nprtestingeb(Result[1,1]',pr);
print "";      
aprtestingeb(Result[1,2:1+rows(matK1)]',pr,matK1);
print ""; 
hprtestingeb(Result[1,2+rows(matK1):1+rows(matK1)+rows(matK2)]',pr,matK2);
print "";
      
apuria=seqa(1989,1/12,rows(z0[1].dataMatrix)); 
graphset;
ylabel("Original Series");          
xlabel("Time");             
if cols(z0[1].dataMatrix)>1;
  _pltype = seqa(1,1,cols(z0[1].dataMatrix));
endif;
xtics(1989,2010,4,1);
_ptek="WL2AR2Series.tkf";
title("Time series of the data");     
xy(apuria,z0[1].dataMatrix);
ret = tkf2eps("WL2AR2Series.tkf", "WL2AR2Series.eps");

apuria=seqa(0,1,rows(z0boots[1].dataMatrix)); 
graphset;
ylabel("Simulated series used in tests");          
xlabel("Time");             
if cols(z0boots[1].dataMatrix)>1;
  _pltype = seqa(1,1,cols(z0boots[1].dataMatrix));
endif;
_ptek="WL2AR2Simuseries.tkf";
title("Time series of the simulated data");     
xy(apuria,z0boots[1].dataMatrix);
ret = tkf2eps("WL2AR2Simuseries.tkf", "WL2AR2Simuseries.eps");

if rows(z0boots[1].dataMatrix)>2000;
  apuria=seqa(0,1,2000); 
  graphset;
  ylabel("Simulated series used in tests");          
  xlabel("Time");             
  if cols(z0boots[1].dataMatrix)>1;
     _pltype = seqa(1,1,cols(z0boots[1].dataMatrix));
  endif;
  _ptek="WL2AR2Simu2000series.tkf";
  title("2000 first obs. of the simulated data");     
  xy(apuria,z0boots[1].dataMatrix[1:2000]);
  ret = tkf2eps("WL2AR2Simu2000series.tkf", "WL2AR2Simu2000series.eps");
endif;

bin=0.10;  
min_z=minc(z0[1].dataMatrix);
max_z=maxc(z0[1].dataMatrix);
Pearson=(seqa(1,1,ceil(1/bin)))/ceil(1/bin)*(max_z-min_z);
Pearson=min_z*ones(rows(Pearson),1)+Pearson;
apuria=seqa(1,1,rows(Pearson)*1000)/1000/rows(Pearson)*(max_z-min_z);
apuria=min_z*ones(rows(apuria),1)+apuria;
apuria=min_z|apuria|max_z+0.0001; 
luvut=counts(z0[1].dataMatrix,Pearson);
luvut=luvut./rows(z0[1].dataMatrix);
graphset;
ylabel("Proportion");  
xlabel("Original series");
xtics(floor(min_z),ceil(max_z),(1/rows(Pearson)),1);
ytics(0,ceil(maxc(100.*luvut))/100,0.1,1);
_pltype={6};
_pcolor={15};
_pxpmax=2;
_pscreen=1;
_ptek="WL2AR2Histo.tkf";
axmargin(1,1,.5,.855);  
title("Histogram of the data");
xy(apuria,0|(luvut).*.ones(1000,1)|0);
ret = tkf2eps("WL2AR2Histo.tkf", "WL2AR2Histo.eps");

apuria=seqa(1989,1/12,rows(mr)); 
graphset;
ylabel("QR Series");          
xlabel("Time");             
if cols(mr)>1;
  _pltype = seqa(1,1,cols(mr));
endif;
xtics(1989,2010,4,1);
_ptek="WL2AR2QRSeries.tkf";
title("Time series of the QRs");     
xy(apuria,mr);
ret = tkf2eps("WL2AR2QRSeries.tkf", "WL2AR2QRSeries.eps");

bin=0.10;  
min_mr=minc(mr);
max_mr=maxc(mr);
Pearson=(seqa(1,1,ceil(1/bin)))/ceil(1/bin)*(max_mr-min_mr);
Pearson=min_mr*ones(rows(Pearson),1)+Pearson;
apuria=seqa(1,1,1000*rows(Pearson))/1000/rows(Pearson)*(max_mr-min_mr);
apuria=min_mr*ones(rows(apuria),1)+apuria;
apuria=min_mr|apuria|max_mr+0.0001; 
luvut=counts(mr,Pearson);
luvut=luvut./rows(mr);
graphset;
ylabel("Proportion");  
xlabel("Quantile Residuals");
xtics(floor(min_mr),ceil(max_mr),(1/rows(Pearson)),1);
ytics(0,ceil(maxc(100.*luvut))/100,0.1,1);
_pltype={6};
_pcolor={15};
_pxpmax=2;
_pscreen=1;
_ptek="WL2AR2QRHisto.tkf";
axmargin(1,1,.5,.855);  
title("Histogram of the Quantile residuals");
xy(apuria,0|(luvut).*.ones(1000,1)|0);
ret = tkf2eps("WL2AR2QRHisto.tkf", "WL2AR2QRHisto.eps");

order_mr=sortc(mr,1);
normals_1=cdfni(seqa(1,1,rows(mr))/(rows(mr)+1));
graphset;
ylabel("Quantile Residuals");  
xlabel("Theoretical Quantiles");
min_val=floor(minc(order_mr[1]|normals_1[1]));
max_val=ceil(maxc(order_mr[rows(mr)]|normals_1[rows(normals_1)]));
line45=seqa(0,1,rows(mr))*(max_val-min_val)/(rows(mr)-1)+min_val;
ytics(min_val,max_val,1,0.5);
xtics(min_val,max_val,1,0.5);
_pltype={2,6};
_pstype={1,0};
_pcolor={15};
_plwidth={20,0};
_pxpmax=0;
_pscreen=1;
_ptek="WL2AR2QR_QQ.tkf";
axmargin(1,1,.5,.855);  
title("Q-Q Plot of Quantile residuals");
xy(normals_1,order_mr~normals_1);
ret = tkf2eps("WL2AR2QR_QQ.tkf","WL2AR2QR_QQ.eps");

bin=1.06*stdc(mr)*rows(mr)^(-1/5);  
R_scale=maxc(abs(mr));
R_domain=seqa(-R_scale,0.01,2*R_scale*100+1);
R_density=ones(rows(R_domain),1);
for i (1,rows(R_domain),1);
  R_density[i]=meanc(pdfn((ones(rows(mr),1)*R_domain[i]-mr)/bin))/bin;
endfor;
R_normal=pdfn((R_domain-meanc(mr))/stdc(mr));
graphset;
ylabel("");  
xlabel("");
xtics(floor(-R_scale),ceil(R_scale),1,1);
ytics(0,ceil(maxc(100.*R_density))/100,0.1,1);
_pltype={6,6};
_pcolor={15,4};
_pxpmax=2;
_plegctl= { 2, 3, 1.7, 4.5 };  /* legend size and locations */
_plegstr= "Kernel estimator \0"\    
          "Scaled normal distribution";
_pscreen=1;
_ptek="WL2AR2QRdensity.tkf";
axmargin(1,1,.5,.855);  
title("N(0,1) Kernel density of the Quantile residuals");
xy(R_domain,R_density~R_normal);
ret = tkf2eps("WL2AR2QRdensity.tkf", "WL2AR2QRdensity.eps");

phi1 = pvUnpack(x0,1);
phi2 = pvUnpack(x0,2);
phi2=phi2|phi1[2:3];
s1 = pvUnpack(x0,3); 
s2 = pvUnpack(x0,4);
beta = pvUnpack(x0,5);
T=rows(z0[1].dataMatrix);
lags=ones(rows(z0[1].dataMatrix),1)~z0[2].dataMatrix~z0[3].dataMatrix;
lags2=ones(rows(z0[1].dataMatrix),1)~z0[3].dataMatrix;
h11=zeros(T,1); h12=h11;
for t (1,T,1);
  h11[t]=(2*pi*det(s1))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi1)'*invpd(s1)*(z0[1].dataMatrix[t]-lags[t,.]*phi1)/2);
  h12[t]=(2*pi*det(s2))^-0.5
       *exp(-(z0[1].dataMatrix[t]-lags[t,.]*phi2)'*invpd(s2)*(z0[1].dataMatrix[t]-lags[t,.]*phi2)/2);
endfor;    
apu=exp(lags2*beta);
p=apu./(ones(rows(apu),1)+apu);
local d,y1,y2,n,muu1,muu2;
n=1;
d=zeros(T,n); y1=zeros(T,1); y2=zeros(T,n); 
d=z0[1].dataMatrix;
muu1=lags*phi1; muu2=lags*phi2;
for t (1,T,1);
  y1[t,1]=ln(g(d[t,.]',muu1[t],s1,muu2[t],s2,p[t]));
endfor;
tau1=p;
tau2=1-p;
for i (1,rows(tau2),1); 
   if tau2[i,1]<=0; tau1[i,1]=tau1[i,1]+tau2[i,1]-1e-8; tau2[i,1]=1e-8; endif; 
endfor;
zhat1=(h11.*tau1)./(h11.*tau1 + h12.*tau2);
zhat2=(1-zhat1);
for i (1,rows(zhat2),1); 
   if zhat2[i,1]<0; zhat1[i,1]=zhat1[i,1]+zhat2[i,1]; zhat2[i,1]=0; endif; 
endfor;
apuria=seqa(1989,1/12,rows(tau1)); 
graphset;
ylabel("Mixing proportions");          
xlabel("Time");             
_pltype ={1,3};
_ptek="WL2AR2MixingSeriesf.tkf";
ytics(0,1,0.1,1);
xtics(1989,2010,3,1);
title("Mixing propostions: Forecast");     
xy(apuria,tau1~(ones(T,1)-tau2));
ret = tkf2eps("WL2AR2MixingSeriesf.tkf", "WL2AR2MixingSeriesf.eps");

m1t=muu1;
m2t=muu2;
conde=zeros(T,1);
condv=conde;
for i (1,T,1);
  conde[i]=tau1[i]*m1t[i]+tau2[i]*m2t[i];
  condv[i]=tau1[i]*s1+tau2[i]*s2
           + tau1[i]*(m1t[i]-conde[i])^2
           + tau2[i]*(m2t[i]-conde[i])^2;
endfor;

graphset;
ylabel("Conditional expectation");          
xlabel("Time");             
_pltype ={1};
_ptek="WL2AR2CondE.tkf";
ytics(0,1,0.1,1);
xtics(1989,2010,3,1);
title("Conditional expectation");     
xy(apuria,conde);
ret = tkf2eps("WL2AR2CondE.tkf", "WL2AR2CondE.eps");

graphset;
ylabel("Conditional variance");          
xlabel("Time");             
_pstype ={1};
_ptek="WL2AR2CondV.tkf";
xtics(1989,2010,3,1);
title("Conditional variance");     
xy(apuria,condv);
ret = tkf2eps("WL2AR2CondV.tkf", "WL2AR2CondV.eps");

graphset;
ylabel("Conditional variance");          
xlabel("Conditional expectation");             
_plctrl = {-1};
_pstype={1};
/*_pltype={0};
_plwidth={2};
_pxpmax=0;
_pscreen=1;*/
_ptek="WL2AR2CondEV.tkf";
title("Observations");     
xy(conde,condv);
ret = tkf2eps("WL2AR2CondEV.tkf", "WL2AR2CondEV.eps");

apuria=seqa(1989,1/12,rows(zhat1)); 
graphset;
ylabel("Mixing proportions");          
xlabel("Time");             
_pltype ={1,3};
_ptek="WL2AR2MixingSeriesi.tkf";
ytics(0,1,0.1,1);
xtics(1989,2010,3,1);
title("Mixing propostions: Inference");     
xy(apuria,zhat1~(ones(T,1)-zhat2));
ret = tkf2eps("WL2AR2MixingSeriesi.tkf", "WL2AR2MixingSeriesi.eps");
 
graphset;
fonts("simplex complex microb simgrma");
ylabel("");          
xlabel("");             
_pcolor={0, 0};
_pltype ={6, 3};
_ptek="WL2AR2All.tkf";
_plegctl= { 3, 5, 300, 400 };  /* legend size and locations */
_plegstr= "y]t[\0"\     /* 4 lines legend text */
          "\204a\201]2,t[*(max(y]t[)-min(y]t[))";
title("");
xtics(1989,2010,3,1);
/*tau1 and zhat1 need to be (0,1)*/ 
min_z=minc(z0[1].dataMatrix);
max_z=maxc(z0[1].dataMatrix);
diff_z=max_z-min_z;    
tau1_scaled=tau1*diff_z+min_z*ones(rows(tau1),1);
zhat1_scaled=zhat1*diff_z+min_z*ones(rows(tau1),1);
xy(apuria,z0[1].dataMatrix~tau1_scaled);
ret = tkf2eps("WL2AR2All.tkf", "WL2AR2All.eps");

KK1=10;
T=rows(mr);
Tboots=rows(mrboots);
k=cols(sboots); covboots=invpd((sboots'*sboots)/Tboots);
GG2=zeros(KK1,k); RR=zeros(KK1,1); guboots=zeros(Tboots-KK1,KK1); help=0;
for t1 (1,(Tboots-KK1),1);
  help=((mrboots[t1+KK1])*mrboots[t1+KK1-1]);
  for i (2,KK1,1);
    help=help~((mrboots[t1+KK1])*mrboots[t1+KK1-i]);
  endfor;
  guboots[t1,.]=help;
endfor;
Sigma=(sboots[KK1+1:Tboots,.]~guboots)'*(sboots[KK1+1:Tboots,.]~guboots)/(Tboots-KK1);
for i1 (1,KK1,1);
    GG2[i1,.]=(meanc((mrboots[1+i1:Tboots]*~mgboots[1:Tboots-i1,.])
              +(mrboots[1:Tboots-i1]*~mgboots[1+i1:Tboots,.])))';
    RR[i1]=meanc(mr[1+i1:T]*~mr[1:T-i1]);    
endfor; 
O2=(GG2*covboots~eye(KK1))*Sigma*(GG2*covboots~eye(KK1))';
autoc=0;
for lk (1,rows(RR),1);
  autoc=autoc|(RR[lk]./sqrt(O2[lk,lk]));
endfor;
autoc=autoc[2:rows(autoc)];
screen on; 
print "ACF of Quantile residuals";
print RR';
print "Scaled ACF of Quantile residuals";
print autoc'; screen off; 
lags=seqa(0,1/100,100*rows(autoc));
autocu=autoc[1]|zeros(99,1);
for kkm (2,rows(autoc),1);
  autocu=autocu|autoc[kkm]|zeros(99,1);
endfor;
autoc=autocu; autoc[1]=0;
ylar=maxc(abs(autoc))+0.02;
critic=0.995;
bounds=(cdfni(critic).*ones(rows(autoc),1))./sqrt(T-KK1);
ylar=maxc(abs(autoc)|bounds)+0.02;
graphset;
  begwind; makewind(9,6.855,0,0,0); makewind(9,6.855,0,0,1); 
  setwind(1);  graphset;    
  ylabel("Autocorrelation");  
  xlabel("Lags");
  title("99% Critical bounds");     
  xtics(0,rows(lags)/100,1,1);
  ytics(-ylar,ylar,(ylar/3),1);
  _pltype={6};
  _pcolor={15};
  _plwidth=15;
  _pxpmax=0;
  _pypmax=2;
  _ptek="WL2AR2AcorQR.tkf";
  axmargin(1,1,.5,.855);  
  xy(lags,autoc);
  nextwind; graphset;    
  axmargin(1,1,.5,.855); 
  ylabel("Autocorrelation");  
  xlabel("Lags"); 
  title("99% Critical bounds");
  xtics(0,rows(lags)/100,1,1);
  ytics(-ylar,ylar,(ylar/3),1);
  _pltype={1,1};
  _pcolor={9,9}; 
  _pxpmax=0;
  _pypmax=2;
  _ptek="WL2AR2AcorQR.tkf";
  xy(lags,bounds~(-bounds));
  endwind;
ret = tkf2eps("WL2AR2AcorQR.tkf", "WL2AR2AcorQR.eps");
  
KK2=10;
T=rows(mr);
Tboots=rows(mrboots);
k=cols(sboots);  covboots=invpd((sboots'*sboots)/Tboots);
GG3=zeros(KK2,k); RRR=zeros(KK2,1); guboots=zeros(Tboots-KK2,KK2);

for i2 (1,KK2,1);
   GG3[i2,.]=2.*meanc(
           (mrboots[KK2+1-i2:Tboots-i2]^2-1)*~mrboots[KK2+1:Tboots]*~mgboots[KK2+1:Tboots,.]+
    ((mrboots[KK2+1:Tboots]^2-1)*~mrboots[KK2+1-i2:Tboots-i2])*~mgboots[KK2+1-i2:Tboots-i2,.])';
   RRR[i2,1]=meanc((mr[KK2+1:T]^2-1)*~(mr[KK2+1-i2:T-i2]^2-1));
   guboots[.,i2]=(mrboots[KK2+1:Tboots]^2-1)*~(mrboots[KK2+1-i2:Tboots-i2]^2-1);   
endfor;
Sigma=(sboots[KK2+1:Tboots,.]~guboots)'*(sboots[KK2+1:Tboots,.]~guboots)/(Tboots-KK2);
O3=(GG3*covboots~eye(KK2))*Sigma*(GG3*covboots~eye(KK2))';                                  

autoc=0;
for lk (1,rows(RRR),1);
  autoc=autoc|(RRR[lk]/sqrt(O3[lk,lk]));
endfor;
autoc=autoc[2:rows(autoc)];
screen on; 
print "ACF of Squared Quantile residuals";
print RRR';
print "Scaled ACF of Squared Quantile residuals";
print autoc'; screen off;
lags=seqa(0,1/100,100*rows(autoc));
autocu=autoc[1]|zeros(99,1);
for kkm (2,rows(autoc),1);
  autocu=autocu|autoc[kkm]|zeros(99,1);
endfor;
autoc=autocu;
autoc[1]=0;
critic=0.995;
bounds=(cdfni(critic).*ones(rows(autoc),1))./sqrt(T-KK2);
ylar=maxc(abs(autoc)|bounds)+0.02;
graphset;
  begwind; makewind(9,6.855,0,0,0); makewind(9,6.855,0,0,1); 
  setwind(1);  graphset;    
  ylabel("Autocorrelation");  
  xlabel("Lags");
  title("99% Critical bounds");
  xtics(0,rows(lags)/100,1,1);
  ytics(-ylar,ylar,(ylar/3),1);
  _pltype={6};
  _pcolor={15};
  _plwidth=15;
  _pxpmax=0;
  _pypmax=2;
  _ptek="WL2AR2AcorQR2.tkf";
  axmargin(1,1,.5,.855);  
  xy(lags,autoc);
  nextwind; graphset;    
  axmargin(1,1,.5,.855); 
  ylabel("Autocorrelation");  
  xlabel("Lags"); 
  title("99% Critical bounds");
  xtics(0,rows(lags)/100,1,1);
  ytics(-ylar,ylar,(ylar/3),1);
  _pltype={1,1};
  _pcolor={9,9}; 
  _pxpmax=0;
  _pypmax=2;
  _ptek="WL2AR2AcorQR2.tkf";
  xy(lags,bounds~(-bounds));
  endwind;
  ret = tkf2eps("WL2AR2AcorQR2.tkf", "WL2AR2AcorQR2.eps");

/* CODES FOR COMPUTING QUANTILE RESIDUAL TESTS */

/* NORMALITY TEST */  
proc ntestingboots(r,rboots,sboots,gboots);    
local Res;    
Res=zeros(1,1);    
Res[1,1]=cdfchic(ntestboots(r,rboots,sboots,gboots),3);
retp(Res);
endp;

proc ntestboots(r,rboots,sboots,gboots);
local G1,O1,T1,T,mv,k,Sigma,Tboots,cov;
T=rows(r);
Tboots=rows(rboots);
k=cols(sboots);
cov=invpd((sboots'*sboots)/Tboots);
G1=(2.*meanc(rboots*~gboots))'|(3.*meanc((rboots^2)*~gboots))'
   |(4.*meanc((rboots^3)*~gboots))';
Sigma=(sboots~(rboots^2-1)~(rboots^3)~(rboots^4-3))'
     *(sboots~(rboots^2-1)~(rboots^3)~(rboots^4-3))/Tboots;
O1=(G1*cov~eye(3))*Sigma*(G1*cov~eye(3))';
mv=(meanc(r^2)-1)|meanc(r^3)|(meanc(r^4)-3);
T1=T*mv'*(invpd(O1))*mv;
retp(T1);
endp;

proc (0) = nprtestingeb(tp,pr); 
print "";
print "*******************************************************************";
print "LM Normality test bootstrapped, r^2-1 included ";
print "  H estim. value   Chi 3 / 3 ";
print tp[1,.]'~pr;   print "";
print "*******************************************************************";
endp;  


/* AUTOCORRELATION TEST */

proc atestingboots(r,rboots,sboots,gboots,matK1); 
local Res,K1,LK1;
K1=1; LK1=rows(matK1); Res=zeros(1,LK1);
for ak1 (1,LK1,1);
   K1=matK1[ak1];
   Res[1,ak1]=cdfchic(atestboots(r,rboots,sboots,gboots,K1),K1);
endfor;
retp(Res);
endp;

proc atestboots(r,rboots,sboots,gboots,K1); 
local O2,G2,T2,T,RR,guboots,k,Sigma,help,Tboots,cov;
T=rows(r);
Tboots=rows(rboots);
k=cols(sboots);
cov=invpd((sboots'*sboots)/Tboots);
G2=zeros(K1,k); RR=zeros(K1,1); guboots=zeros(Tboots-K1,K1); help=0;
for t1 (1,(Tboots-K1),1);
  help=((rboots[t1+K1])*rboots[t1+K1-1]);
  for i (2,K1,1);
    help=help~((rboots[t1+K1])*rboots[t1+K1-i]);
  endfor;
  guboots[t1,.]=help;
endfor;
Sigma=(sboots[K1+1:Tboots,.]~guboots)'*(sboots[K1+1:Tboots,.]~guboots)/(Tboots-K1);
for i1 (1,K1,1);
    G2[i1,.]=(meanc((rboots[1+i1:Tboots]*~gboots[1:Tboots-i1,.])
              +(rboots[1:Tboots-i1]*~gboots[1+i1:Tboots,.])))';          
    RR[i1]=meanc(r[1+i1:T]*~r[1:T-i1]);    
endfor; 
O2=(G2*cov~eye(K1))*Sigma*(G2*cov~eye(K1))';
T2=(T-K1)*(RR')*(invpd(O2))*RR;
retp(T2);
endp;

proc (0) = aprtestingeb(tp,pr,matK1);
local LK1,label;
LK1=rows(matK1);
print "";
print "*******************************************************************";
print "LM AUTOCORRELATION TESTS with bootstrapped omega";
label="LM Autocorrelation tests";
call pTEST(label,matK1,tp[1:LK1,.],pr);
print "";
print "*******************************************************************";
endp;

/* CONDITIONAL HETEROSKEDASTICITY TEST */

proc htestingboots(r,rboots,sboots,gboots,matK2); 
local Res,K2,LK2;
K2=1; LK2=rows(matK2); Res=zeros(1,LK2);
for ak2 (1,LK2,1);
   K2=matK2[ak2];
   Res[1,ak2]=cdfchic(htestboots(r,rboots,sboots,gboots,K2),K2);
endfor;
retp(Res);
endp;

proc htestboots(r,rboots,sboots,gboots,K2); 
local O3,G3,T3,T,RRR,k,Sigma,guboots,Tboots,cov;
T=rows(r); Tboots=rows(rboots);
k=cols(sboots); cov=invpd((sboots'*sboots)/Tboots);
G3=zeros(K2,k); RRR=zeros(K2,1);
guboots=zeros(Tboots-K2,K2);
for i2 (1,K2,1);
   G3[i2,.]=2.*meanc(
           (rboots[K2+1-i2:Tboots-i2]^2-1)*~rboots[K2+1:Tboots]*~gboots[K2+1:Tboots,.]+
    ((rboots[K2+1:Tboots]^2-1)*~rboots[K2+1-i2:Tboots-i2])*~gboots[K2+1-i2:Tboots-i2,.])';
   RRR[i2,1]=meanc((r[K2+1:T]^2-1)*~(r[K2+1-i2:T-i2]^2-1));
   guboots[.,i2]=(rboots[K2+1:Tboots]^2-1)*~(rboots[K2+1-i2:Tboots-i2]^2-1);   
endfor;
Sigma=(sboots[K2+1:Tboots,.]~guboots)'*(sboots[K2+1:Tboots,.]~guboots)/(Tboots-K2);
O3=(G3*cov~eye(K2))*Sigma*(G3*cov~eye(K2))';                                  
T3=T*(RRR')*(invpd(O3))*RRR;
retp(T3);
endp;

proc (0) = hprtestingeb(tp,pr,matK2);
local LK2,label; LK2=rows(matK2);
print "";
print "*******************************************************************";
print "LM HETEROSCEDASTICITY TEST with bootstrapped omega";
label="LM Heteroscedasticity test";
call pTEST(label,matK2,tp[1:LK2,.],pr);
print "";
print "*******************************************************************";
endp;


proc (0) = pTEST(label,matK,tp,pr);
local LK;
LK=rows(matK);
for i(1,rows(matK),3);  
  print label; 
  if LK>=3; 
    print "          Lags " ftos(matK[i],"%*.*lf",1,0) "           Lags " 
     ftos(matK[i+1],"%*.*lf",1,0) "           Lags " ftos(matK[i+2],"%*.*lf",1,0) 
     "         Chi lags";
    print tp[i:i+2,.]'~pr;
    print "";    
  endif;
  if LK==2;
    print "          Lags " ftos(matK[i],"%*.*lf",1,0) "           Lags " 
     ftos(matK[i+1],"%*.*lf",1,0) "        Chi lags"; 
    print tp[i:i+1,.]'~pr;
    print ""; 
  endif;
  if LK==1;
    print "          Lags " ftos(matK[i],"%*.*lf",1,0)  "        Chi lags";
    print tp[i,.]'~pr;
    print ""; 
  endif;
  LK=LK-3;
endfor; 
endp; 

proc resid(&g,&cumulg,data,mu,sigma,nu,omega,p);
local g:proc, cumulg:proc;
local eps;
eps=10^-10;
if rows(data)<=0;
  retp (error(0));
elseif rows(data)==1;
  if cumulg(data,mu,sigma,nu,omega,p)>=1;
    retp (1-eps);
  else;
    if cumulg(data,mu,sigma,nu,omega,p)<=0;
      retp (0+eps);
    else;
      retp (cdfni(cumulg(data,mu,sigma,nu,omega,p)));
    endif;
  endif; 
else;
  local x1,x2,mu1,mu2,sigma1,sigma2,nu1,nu2,omega1,omega2,p2,e;
  x1=data[1,1];
  x2=data[2:rows(data),1];
  mu1=mu[1,1];
  nu1=nu[1,1];
  sigma1=sigma[1,1];
  omega1=omega[1,1];
  mu2=mu[2:rows(mu),1]+sigma[2:rows(sigma),1]*(x1-mu1)/sigma1;
  nu2=nu[2:rows(nu),1]+omega[2:rows(omega),1]*(x1-nu1)/omega1;
  sigma2=sigma[2:rows(sigma),2:cols(sigma)]
         -sigma[2:rows(sigma),1]*sigma[1,2:cols(sigma)]/sigma1;
  omega2=omega[2:rows(omega),2:cols(omega)]
         -omega[2:rows(omega),1]*omega[1,2:cols(omega)]/omega1;
  p2=p*g(x1,mu1,sigma1,nu1,omega1,1)/g(x1,mu1,sigma1,nu1,omega1,p);  
  e=zeros(rows(data),1);
  eps=10^-10;
  if cumulg(x1,mu1,sigma1,nu1,omega1,p)>=1;
    e[1,1]=(1-eps);
  else;
    if cumulg(x1,mu1,sigma1,nu1,omega1,p)<=0;
      e[1,1]=(0+eps);
    else;
      e[1,1]=cdfni(cumulg(x1,mu1,sigma1,nu1,omega1,p));
    endif;
  endif;
  e[2:rows(data),1]=resid(&g,&cumulg,x2,mu2,sigma2,nu2,omega2,p2);
  retp (e);
endif;
endp;

proc g(data,m1,s1,m2,s2,p);
local xx1,xx2, k, y;
xx1=(data-m1)'*invpd(s1)*(data-m1);
xx2=(data-m2)'*invpd(s2)*(data-m2);
k=rows(data); 
y=    p*det(s1)^(-0.5)*(2*pi)^(-k/2)*exp(-0.5*xx1)
 +(1-p)*det(s2)^(-0.5)*(2*pi)^(-k/2)*exp(-0.5*xx2);
retp (y);
endp;
 
proc cumulg(data,m1,s1,m2,s2,p);
local x1,x2;
x1=(data-m1)/sqrt(s1);
x2=(data-m2)/sqrt(s2); 
retp (p*cdfn(x1)+(1-p)*cdfn(x2));
endp;