/*    22.05.2011  at 16:00
**
**    COMPUTES various test statistics:
**    Normality tests
**    Autocorrelation tests
**    Heteroscedasticity tests
**
**    With large sample based estimators for covariance matrix
**    in the test statistics
**
**    PRINTS the tests commands 
**
*/

/*
*****************************************************************************
**  UNIVARIATE SERIES - TESTS                                                       ** 
*****************************************************************************
*/

/* Normality test - TYPE 1 - (r^2-1) included
**
**   FORMAT    T1 = ntesttype1(r,rboots,sboots,gboots);     
**               
**   INPUT      r - the vector of quantile residuals  (T x 1)  
**         rboots - simulated vector of QRs using estimated model  (Tboots x 1)        
**         sboots - similated matrix of scores of QRs using estimated model 
**                 (Tboots x k)                                  
**         gboots - simulated matrix of derivatives of QRs using estimated models 
**                 (Tboots x k)                                   
**            
** 
**   OUTPUT    T1 - scalar value of the test statistic
**
*/          

proc ntesttype1(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;
/****************************************************************************/

/* Normality test - TYPE 2
**
**   FORMAT    T1 = ntesttype2(r,rboots,sboots,gboots);     
**               
**   INPUT      r - the vector of quantile residuals  (T x 1)  
**         rboots - simulated vector of QRs using estimated model  (Tboots x 1)        
**         sboots - similated matrix of scores of QRs using estimated model 
**                 (Tboots x k)                                  
**         gboots - simulated matrix of derivatives of QRs using estimated models 
**                 (Tboots x k)                                   
**     
** 
**   OUTPUT    T1 - scalar value of the test statistic
**
*/          

proc ntesttype2(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=(3.*meanc((rboots^2)*~gboots))'
   |(4.*meanc((rboots^3)*~gboots))';
Sigma=(sboots~(rboots^3)~(rboots^4-3))'
     *(sboots~(rboots^3)~(rboots^4-3))/Tboots;
O1=(G1*cov~eye(2))*Sigma*(G1*cov~eye(2))';
mv=meanc(r^3)|(meanc(r^4)-3);
T1=T*mv'*(invpd(O1))*mv;
retp(T1);
endp;
/****************************************************************************/


/* Autocorrelation test 
**
**   FORMAT    T2 = atest(r,rboots,sboots,gboots,K1);     
**               
**   INPUT      
**              r - the vector of quantile residuals  (T x 1)        
**         rboots - simulated vector of QRs using estimated model  (Tboots x 1)        
**         sboots - similated matrix of scores of QRs using estimated model 
**                 (Tboots x k)                                  
**         gboots - simulated matrix of derivatives of QRs using estimated models 
**                 (Tboots x k)                                    
**             K1 - the number of lags used in the test 
**           
** 
**   OUTPUT    T2 - scalar value of the test statistic
**
*/    

proc atest(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;

/****************************************************************************/

/****************************************************************************/

/* Heteroscedasticity test 
**
**   FORMAT    T3 = htest(r,rboots,sboots,gboots,K2);     
**               
**   INPUT     
**              r - the vector of quantile residuals  (T x 1)        
**         rboots - simulated vector of QRs using estimated model  (Tboots x 1)        
**         sboots - similated matrix of scores of QRs using estimated model 
**                 (Tboots x k)                                  
**         gboots - simulated matrix of derivatives of QRs using estimated models 
**                 (Tboots x k)    
**             K2 - the number of lags used in the test 
**                     
**   OUTPUT    T3 - scalar value of the test statistic
**
*/    

proc htest(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;
/****************************************************************************/


/*
*****************************************************************************
**  MULTIVARIATE SERIES - TESTS                                            **
*****************************************************************************
*/


/* Multinormality test - TYPE 1 - (r^2-1) included
**
**   FORMAT    T1 = multintesttype1(r,rboots,sboots,gboots);     
**               
**   INPUT      r - the matrix of multivariate quantile residuals  (T x n)        
**         rboots - simulated matrix of multivariate quantile residuals 
**                  using estimated model (Tboots x k)    
**         sboots - simulated matrix of scores of multivariate quantile residuals 
**                  using estimated model (Tboots x k)                                   
**         gboots - simulated matrix of derivatives of multivariate quantile residuals
**                  (Tboots x nk) 
** 
**   OUTPUT    T1 - scalar value of the test statistic
**
*/      

proc multintesttype1(r,rboots,sboots,gboots); 
local G1,O1,T1,T,n,gu,raboots,k,Sigma,Tboots,cov;
T=rows(r);
k=cols(sboots);
Tboots=rows(rboots);
n=cols(r);
cov=invpd((sboots'*sboots)/Tboots);
/* Martix G1 has to be n4*k                                                 */
/* meanc(g[((i-1)*T+1):i*T,.])' is (1 x k matrix)                           */
G1=2.*meanc(rboots[.,1]*~gboots[.,1:k])'|
  (3.*meanc((rboots[.,1]^2)*~gboots[.,1:k]))'|
  (4.*meanc((rboots[.,1]^3)*~gboots[.,1:k]))';
/* Vector gu has to be n4*1                                                 */
gu=meanc(r[.,1]^2-1)|meanc(r[.,1]^3)|meanc(r[.,1]^4-3);
raboots=(rboots[.,1]^2-1)~(rboots[.,1]^3)~(rboots[.,1]^4-3);
for i (2,n,1);
  G1=G1|
     (2.*meanc(rboots[.,i]*~gboots[.,((i-1)*k+1):i*k]))'|
     (3.*meanc((rboots[.,i]^2)*~gboots[.,((i-1)*k+1):i*k]))'|
     (4.*meanc((rboots[.,i]^3)*~gboots[.,((i-1)*k+1):i*k]))';
  gu=gu|
     meanc(r[.,i]^2-1)|meanc(r[.,i]^3)|meanc(r[.,i]^4-3);
  raboots=raboots~(rboots[.,i]^2-1)~(rboots[.,i]^3)~(rboots[.,i]^4-3);
endfor;
Sigma=(sboots~raboots)'*(sboots~raboots)/Tboots;
O1=(G1*cov~eye(3*n))*Sigma*(G1*cov~eye(3*n))';
T1=T*gu'*(invpd(O1))*gu;
retp(T1);
endp;
/****************************************************************************/


/* Multinormality test - TYPE 2 
**
**   FORMAT    T1 = multintesttype2(r,rboots,sboots,gboots);     
**               
**   INPUT      r - the matrix of multivariate quantile residuals  (T x n)        
**         rboots - simulated matrix of multivariate quantile residuals 
**                  using estimated model (Tboots x n)    
**         sboots - simulated matrix of scores of multivariate quantile residuals 
**                  using estimated model (Tboots x k)                                   
**         gboots - simulated matrix of derivatives of multivariate quantile residuals
**                  (Tboots x nk) 
** 
**   OUTPUT    T1 - scalar value of the test statistic
**
*/      

proc multintesttype2(r,rboots,sboots,gboots); 
local G1,O1,T1,T,n,gu,raboots,k,Sigma,Tboots,cov;
T=rows(r);
k=cols(sboots);
Tboots=rows(rboots);
n=cols(r);
cov=invpd((sboots'*sboots)/Tboots);
/* Martix G1 has to be n2*k                                                 */
/* meanc(g[((i-1)*T+1):i*T,.])' is (1 x k matrix)                           */
G1=(3.*meanc((rboots[.,1]^2)*~gboots[.,1:k]))'|
   (4.*meanc((rboots[.,1]^3)*~gboots[.,1:k]))';
/* Vector gu has to be n4*1                                                 */
gu=meanc(r[.,1]^3)|meanc(r[.,1]^4-3);
raboots=(rboots[.,1]^3)~(rboots[.,1]^4-3);
for i (2,n,1);
  G1=G1|
     (3.*meanc((rboots[.,i]^2)*~gboots[.,((i-1)*k+1):i*k]))'|
     (4.*meanc((rboots[.,i]^3)*~gboots[.,((i-1)*k+1):i*k]))';
  gu=gu|
     meanc(r[.,i]^3)|meanc(r[.,i]^4-3);
  raboots=raboots~(rboots[.,i]^3)~(rboots[.,i]^4-3);
endfor;
Sigma=(sboots~raboots)'*(sboots~raboots)/Tboots;
O1=(G1*cov~eye(2*n))*Sigma*(G1*cov~eye(2*n))';
T1=T*gu'*(invpd(O1))*gu;
retp(T1);
endp;
/****************************************************************************/


/* Multivariate Autocorrelation test 
**
**   FORMAT    T2 = multiatest(r,rboots,sboots,gboots,K1);     
**               
**   INPUT      
**              r - the matrix of multivariate quantile residuals  (T x n)
**         rboots - simulated matrix of multivariate quantile residuals 
**                  using estimated model (Tboots x n)    
**         sboots - simulated matrix of scores of multivariate quantile residuals 
**                  using estimated model (Tboots x k)                                   
**         gboots - simulated matrix of derivatives of multivariate quantile residuals
**                  (Tboots x nk) 
**             K1 - the number of lags used in the test 
**           
**   OUTPUT    T2 - scalar value of the test statistic
**
*/    

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

/* Multivariate Heteroscedasticity test
**
**   FORMAT    T3 = multihtest(r,rboots,sboots,gboots,K2);     
**               
**   INPUT      
**              r - the vector of multivariate quantile residuals  (T x n)
**         rboots - simulated matrix of multivariate quantile residuals 
**                  using estimated model (Tboots x n)    
**         sboots - simulated matrix of scores of multivariate quantile residuals 
**                  using estimated model (Tboots x k)                                   
**         gboots - simulated matrix of derivatives of multivariate quantile residuals
**                  (Tboots x nk) 
**             K2 - the number of lags used in the test 
**           
** 
**   OUTPUT    T3 - scalar value of the test statistic
**
*/    

proc multihtest(r,rboots,sboots,gboots,K2); 
local O3,G3,T3,T,n,k,gu,guboots,RRR,help,helpboots,Sigma,Tboots,cov;
T=rows(r);
Tboots=rows(rboots);
n=cols(r);
k=cols(sboots);
cov=invpd((sboots'*sboots)/Tboots);
gu=zeros(T-K2,K2*n^2);
help=0;
for t1 (1,(T-K2),1);
  help=((r[t1+K2,.]^2-1)'*(r[t1+K2-1,.]^2-1));
  for i (2,K2,1);
    help=help~((r[t1+K2,.]^2-1)'*(r[t1+K2-i,.]^2-1));
  endfor;
  gu[t1,.]=vecr(help')'; 
endfor;
RRR=meanc(gu);
guboots=zeros(Tboots-K2,K2*n^2);
helpboots=0;
for t1 (1,(Tboots-K2),1);
  helpboots=((rboots[t1+K2,.]^2-1)'*(rboots[t1+K2-1,.]^2-1));
  for i (2,K2,1);
    helpboots=helpboots~((rboots[t1+K2,.]^2-1)'*(rboots[t1+K2-i,.]^2-1));
  endfor;
  guboots[t1,.]=vecr(helpboots')'; 
endfor;
G3=zeros(K2*n^2,k);
for j (1,K2,1);
  for b (1,n,1);
    for a (1,n,1);
      G3[((j-1)*n^2+(b-1)*n+a),.]=
        2.*meanc(
           (rboots[(K2+1-j):(Tboots-j),b]^2-1)*~rboots[(K2+1):Tboots,a]
           *~gboots[(K2+1):Tboots,((a-1)*k+1):a*k]+
           (rboots[K2+1:Tboots,a]^2-1)*~rboots[(K2+1-j):(Tboots-j),b]
           *~gboots[(K2+1-j):(Tboots-j),((b-1)*k+1):b*k]
        )';   
    endfor;
  endfor;
endfor;
Sigma=(sboots[K2+1:Tboots,.]~guboots)'*(sboots[K2+1:Tboots,.]~guboots)/(Tboots-K2);
O3=(G3*cov~eye(K2*n^2))*Sigma*(G3*cov~eye(K2*n^2))';                                  
T3=T*(RRR')*(invpd(O3))*RRR;
retp(T3);
endp;

/****************************************************************************/


/*   Prints test statistics p-values 
**   3 tests at the same row 
**  
**   FORMAT    pTEST(matK,tp);     
**               
**   INPUT       matK - vector, containing the lag values used in calculation
**                      of the test statistics
**                 tp - matrix ( rows(matK) x 1 ) contains the estimated 
**                      p-values of the test statistics
**   OUTPUT    returns no value
**
*/            

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

    print tp[i,.]';
    print ""; 
  endif;
  LK=LK-3;
endfor; 
  
endp;      
/*******************************************************************************/
/*******************************************************************************/

/*
**   FORMAT    prtesting(r,rboots,sboots,gboots,matK1,matK2);     
**               
**   INPUT       see proc's "ntesttype1", "ntesttype2",
**               "atest" & "htest".
**
**   OUTPUT    returns no value
*/

proc (0) = prtesting(r,rboots,sboots,gboots,matK1,matK2); 
print "";
print "*******************************************************************";
print "p-value of the LM Normality test TYPE 1 (r^2-1 included)";
print cdfchic(ntesttype1(r,rboots,sboots,gboots),3);   
print "";
print "p-value of the LM Normality test TYPE 2";
print cdfchic(ntesttype2(r,rboots,sboots,gboots),2);   
print "";
print "*******************************************************************";
local Res1,K1,LK1;
K1=1;
LK1=rows(matK1);
Res1=zeros(1,LK1);
for ak1 (1,LK1,1);
   K1=matK1[ak1];
   Res1[1,ak1]=cdfchic(atest(r,rboots,sboots,gboots,K1),K1);
endfor;
print "";
print "*******************************************************************";
print "LM AUTOCORRELATION TESTS";
call pTEST(matK1,Res1');
print "";
print "*******************************************************************";
local Res2,K2,LK2;
K2=1;
LK2=rows(matK2);
Res2=zeros(1,LK2);
for ak2 (1,LK2,1);
   K2=matK2[ak2];
   Res2[1,ak2]=cdfchic(htest(r,rboots,sboots,gboots,K2),K2);
endfor;
print "";
print "*******************************************************************";
print "LM HETEROSCEDASTICITY TEST ";
call pTEST(matK2,Res2');
print "";
print "*******************************************************************";
endp;
/*******************************************************************************/

/*
**      MULTIVARIATE TESTS
*/

/*******************************************************************************/

/*
**   FORMAT    multiprtesting(r,rboots,sboots,gboots,matK1,matK2);     
**               
**   INPUT     See proc's "multintesttype1", "multintesttype2",
**             "multiatest" & "multihtest"   
** 
**   OUTPUT    returns no value
**
*/        
proc (0) = multiprtesting(r,rboots,sboots,gboots,matK1,matK2);
local n;
n=cols(r);    
print "*******************************************************************";
print "";
print "Multinormality test TYPE 1 (r^2-1 included) ";
print cdfchic(multintesttype1(r,rboots,sboots,gboots),3*n);
print "";
print "Multinormality test TYPE 2 ";
print cdfchic(multintesttype2(r,rboots,sboots,gboots),2*n); 
print "";
print "*******************************************************************";
local Res1,K1,LK1;
K1=1;
LK1=rows(matK1);
Res1=zeros(1,LK1);
for ak1 (1,LK1,1);
  K1=matK1[ak1];      
  Res1[1,ak1]=cdfchic(multiatest(r,rboots,sboots,gboots,K1),K1*n^2);
endfor;
print "*******************************************************************";
print ""; 
print "LM MULTIVARIATE AUTOCORRELATION TESTS";
call pTEST(matK1,Res1');
print "";
print "*******************************************************************";
print "";
local Res2,K2,LK2;
K2=1;
LK2=rows(matK2);
Res2=zeros(1,LK2);
for ak2 (1,LK2,1);
   K2=matK2[ak2];     
   Res2[1,ak2]=cdfchic(multihtest(r,rboots,sboots,gboots,K2),K2*n^2);
endfor;
print "*******************************************************************";
print "";
print "LM MULTIVARIATE HETEROSCEDASTICITY TESTS";
call pTEST(matK2,Res2'); print "";    
print "*******************************************************************";
print "";
endp;