/*    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;
/****************************************************************************/


/*   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;
/*******************************************************************************/