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

	Some general purpose procedures that are used by MLAGNR.G
	in the computations.
	Created by Antti Ripatti, 26.10.1998

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

proc (1)=min(x);
     local p,i,j,r;

     x=vec(x);
     i=1; j=1; p=x[1];
     r=rows(x);
     do while i<=r;
        do while j<=r;
            if x[j]<=p;
                p=x[j];
            endif;
            j=j+1;
        endo;
        i=i+1;
     endo;
     retp(p);
endp;


proc (1)=max(x);
     local p,i,j,r;

     x=vec(x);
     i=1; j=1; p=x[1];
     r=rows(x);
     do while i<=r;
        do while j<=r;
            if x[j]>=p;
                p=x[j];
            endif;
            j=j+1;
        endo;
        i=i+1;
     endo;
     retp(p);
endp;

/* ACF
**
** Purpose:   compute autocorrelation function
**
** Format:    x=ACF(y,maxl);
**
** Input:     y     time series vector
**            maxl  scalar, max lag length
**
** Output:    x     vector, autocorrelations
**
** Created:   11.4.94 by Antti Ripatti
** Revised:
**
*/
proc (1)=acf(y,maxl);
     local t,v,yt,ytk,vk,p,i;

     y=y-meanc(y);
     t=rows(y);
     v=y'y/t;
     p=zeros(maxl,1);
     i=1;
     do while i <= maxl;
        yt=y[1:t-i,1];
        ytk=y[1+i:t,1];
        vk=yt'ytk/t;
        p[i,1]=vk/v;
        i=i+1;
     endo;
     retp(p);
endp;


/* ADJOBS
**
** Purpose:   returns the last nobs rows of the matrix
**
** Format:    y = ADJOBS(x,nobs);
**
** Input:     x     T x K matrix, time serie
**            nobs  scalar, the desired number of rows of x
**
** Output:    y     T-nobs x K matrix, where the first T-nobs rows are deleted
**            if nobs>T then prints error message and returns scalar -1 
**
** Created:   15.4.96 by AR
** Revised:
**
*/
proc (1) = AdjObs(x,nobs);
    local xrows, retvalue;

    xrows = rows(x);
    if xrows < nobs;
        errorlog "Error in calling procedure ADJOBS. Number of rows desired";
        errorlog "exceeds the number of rows available. Function returns -1";
        retvalue = -1;
    else;
        retvalue = x[xrows-nobs+1:xrows,.];
    endif;
    retp(retvalue);
endp;


proc (1) = ao(x,nobs);

    retp(adjobs(x,nobs));
endp;

/* DIFF */
/* Calculates first differeces
** >> Does NOT add a miss. value in place of the first observation
*/
FN diff(X) = trimr(X,1,0)-trimr(X,0,1);



proc (0)=dims(x);
    print ": " rows(x) " x " cols(x);
endp;


/* DEMEAN
**
** Purpose:   to deseasonalize, to demean the original series
**
** Format:    y = demean(X,S);
**
** Input:     X     T x k   matrix of matrix to be demeaned
**            S     T x s   matrix of deterministic components
**            
** Output:    y     T x k   matrix of demeaned X
**
** Created:   16.9.1995 by Antti Ripatti
** Revised:
**
*/
proc (1)=demean(X,S);
    local y, nobs;

    X=packr(X);
    S=packr(S);
    nobs = min(rows(X)~rows(S));
    S=ao(S,nobs);
    X=ao(X,nobs);
    y =  X - S*inv(S'S)*S'X;
    retp(y);
endp;

/* DESEAS
**
** Purpose:   to deseasonalize, to demean the original series
**
** Format:    y = deseas(X,S);
**
** Input:     X     T x k   matrix of matrix to be demeaned
**            S     T x s   matrix of deterministic components
**            
** Output:    y     T x k   matrix of demeaned X
**
** Created:   16.9.1995 by Antti Ripatti
** Revised:
**
*/
proc (1)=DeSeas(X,S);
    local y, nobs, beta;

    X=packr(X);
    S=packr(S);
    nobs = min(rows(X)~rows(S));
    S=ao(S,nobs);
	S=ones(nobs,1)~S;
    X=ao(X,nobs);
	beta = inv(S'S)*S'X;
    y =  X - S*beta + ones(nobs,1)*beta[1,.];
    retp(y);
endp;

proc (1)=createZ0(Z,k);
     local T, dZ;

     T=rows(Z);
     dZ=diff(Z);
     T=T-1; /* note that rows(dZ) = rows(Z)-1; */
     retp(dZ[k:T,.]);
endp;

proc (1)=createZ1(Z,k,D);
     local T, i, dZ, Z1, tmp;

     T=rows(Z);
     dZ=diff(Z);
     T=T-1;
     Z1=dZ[k-1:T-1,.];
     i=2;
     do while i<=(k-1);
        tmp=Z1~dZ[k-i:T-i,.];
        Z1=tmp;
        i=i+1;
     endo;
     if rows(D)>1;
        Z1=Z1~D[k:T,.];
     endif;
     retp(Z1);
endp;

proc (1)=createZk(Z,k,const);
     local T;

     T=rows(Z)-1;
     if const==0;
        retp(Z[k:T,.]);
     else;
        retp(Z[k:T,.]~ones(T-(k-1),1));
     endif;
endp;


proc (4)=urrr(Z0,Z1,Zk,p);
     local T, i, R0, Rk, S00, S0k, Skk, Sh, SS, l, V, beta, alfa, g,
           pii, Lr, r, tracet;

     T=rows(Z0);
     R0=Z0 - Z1*inv(Z1'Z1)*Z1'Z0;
     Rk=Zk - Z1*inv(Z1'Z1)*Z1'Zk;
     S00=moment(R0,0)/T;
/* "-Tln(|S00|)";;-T*ln(det(S00)); */
     S0k=(R0'Rk)/T;
     Skk=moment(Rk,0)/T;
     clear R0, Rk;
     Sh=inv(chol(Skk));
/*
     matrix to be solved to achieve ML estimator
*/
     SS=Sh'S0k'invpd(S00)*S0k*Sh;
     {l, V}=eighv(SS);
     l=rev(l);
     V=(rev(V'))';
     beta=Sh*V;
     alfa=S0k*beta;
     pii=alfa*beta';
     beta=beta./beta[1,.];
     alfa=(pii'/beta)';
     l=l[1:p];
     beta=beta[.,1:p];
     alfa=alfa[1:p,1:p];
     g=rows(l);
     Lr=-T*ln(ones(g,1)-l);
     r=0; tracet=zeros(p,1);
     do while r<p;
        tracet[r+1]=sumc(Lr[r+1:p]);
        r=r+1;
     endo;
     retp(l,beta,alfa,tracet);
endp;


/* TS
**
** Purpose:    graph time series picture
**
** Format:     TS(x,b,p);
**
** Input:     x  time series vector
**            b  scalar, starting year of x
**            p  scalar, periodicity of x
**
** Output:    none
**
** Created:  11.4.94 by Antti Ripatti
** Revised
**
*/
proc (0)=ts(x,b,p);
    local rx, toler;

    toler=1E-300;
    rx=rows(x);
    _pdate="";
    _pnum={1 2};
    _pypmax=3;
    _pysci=3;
/*    _plwidth=2; */
    xtics(b,b+ceil(rx/p),1,p);
    if abs(min(x)-max(x)) < toler;
        ytics(min(x)-toler,max(x)+toler,toler,3);
    else;
        ytics(min(x),max(x),(max(x)-min(x))/5.5,5);
    endif;
    xy(seqa(b,1/p,rx),x);
endp;

/* create r x c matrix of missing values */
proc (1) = crmiss(r,c);
    retp(miss(ones(r,c),ones(r,c)));
endp;




/* DISPLAYRESULTS
**
** Purpose: to display and format nicely parameter names, values, std.errs,
**          t-values and p-values
**
** Format:  DisplayResults(Covariance, Params, ParameterNames, NObs);
**
** Input:   Covariance  p x p matrix covariances
**          Params      p x 1 matrix of parameter values
**          ParameterNames p x 1 matrix of strings, parameter names
**          NObs        scalar, number of observations
**
** Output:  print out
**
** Created: 25.9.1996 by Antti Ripatti
** Revised:
**
*/
proc (0)=DisplayResults(Covariance, Params, ParameterNames, NObs);
    local SDs, mask, fmt, PrtMat;
    
    SDs     = sqrt(abs(diag(Covariance))); /* note that using ABS is ILLEGAL, variance should be positive */
    if sumc(packr(diag(Covariance)).<0) > 0;
        print "WARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWARNINGWA";
        print " Some of the diagonal elements of the covariance-variance matrix";
        print "  are NEGATIVE! Absolute values used in computing std. errors! ";
        print "   **********************************************************  ";
    endif;
    mask    = zeros(rows(Params),1)~ones(rows(Params),4);
    mask    = 0~1~1~1~1;
    let fmt[5,3]  = "-*.*s" 8 8 "*.*lf" 15 6 "*.*lf" 15 6 "*.*lf" 15 6 "*.*lf" 15 6;
    print "Parameter         value       std.err.    t-statistic        p-value";
    print "--------------------------------------------------------------------";
    prtmat = ParameterNames~Params~SDs~(Params./SDs)~(2*cdftc(abs(Params./SDs),nobs));
    call printfm(prtmat,mask,fmt);
    print;
endp;



proc (7)=resdiag(res, T, ssors, lags);
     local sk, ek, jb, hq, q1, q4, q12;

     sk=skewness(res);
     ek=kurtosis(res);
     q1=T*(T+2)*acf(res,1)^2/(T-1);
     q4=T*(T+2)*sumc(((acf(res,4))^2)/(T-1));
     q12=T*(T+2)*sumc(((acf(res,12))^2)/(T-1));
     jb=jbnorm(sk,ek,T,ssors);
     hq=ln(sumc(res^2)/(T-1))+lags*2*ln(ln(T))/T;
     retp(sk,ek,jb,hq,q1,q4,q12);
endp;


/* DResDiag
**
** Purpose:   Compute and display some residual statistics
**
** Format:    DResDiag(res,nobs, Ssors, lags, CLab);
**
** Input:     res   T x n matrix of residuals
**            nobs  scalar, number of observations
**            Ssors scalar, number of regressors
**            lags  scalar, number of lags in the model
**            Clab  n x 1  array of strings
**
** Output:    returns 1 when no problems, 0 otherwise
**
** Created:   6.11.1996 by Antti Ripatti
** Revised:
**
*/
proc (1)=DResDiag(res, NObs, ssors, lags, clab);
     local sk, ek, jb, hq, q1, q4, q12, i, head, cvm, stdm, Eqs, tmp;

     i=1;
     Eqs = cols(res);
     do while i<=Eqs;
        {sk,ek,jb,hq,q1,q4,q12}=resdiag(res[.,i],Nobs,ssors, lags);
        print ;
        format /m0 /rd 10,5;
        print "EQUATION: " $CLab[i];
        print "    Hannan & Quinn information criteria  :" hq;
        print "    Mean                                 :" meanc(res[.,i]);
        tmp = res[.,i] - meanc(res[.,i]);
        print "    Variance                             :" moment(tmp,0)/nobs;
        print "    Skewness                             :" sk;
        print "    Excess Kurtosis                      :" ek;
        print "    Normality                            :" jb "     (" cdfchic(jb,2) ")";
        print "    Box-Pierce Q 1  lag                  :" q1 "     (" cdfchic(q1,1) ")";
        print "               Q 4  lags                 :" q4 "     (" cdfchic(q4,4) ")";
        print "               Q 12 lags                 :" q12 "     (" cdfchic(q12,12) ")";
        print "    Residual autocorrelation             :";;
        format /m1 /rds 5,2;
        (seqa(1,1,12)~acf(res[.,i],12))';
        format /m0 /rd 10,5;
        print " ";print " ";
        i=i+1;
     endo;
     print "Residual variance - correlation matrix:";
     cvm=corrx(res);
     stdm=diag(moment(res,0))-ones(cols(res),1); 
     format /m1 /rd 10,5;
     print cvm+diagrv(zeros(rows(stdm),rows(stdm)),stdm);
     print;
     retp(1);
endp;


/* RESGRAPH
**
** Purpose: Is to graph residuals and the residual autocorrelation function   
**
** Format:      call resgraph(res,year,period,periodicity,lags,tofile);   
**
** Input:       res         T x 1 vector of residuals,
**              year        scalar of starting year
**              period      scalar of startin period
**              periodicity scalar of periodicity of the data
**              lags        scalar, number of autocorrelation terms
**              tofile      string, 0 if displayed on screen, otherwise the name of EPS-file
**
** Output:      2 by 1 graph
**
** Created:     25.4.1996 by Antti Ripatti
** Revised:
**
*/
proc (0)=resgraph(res,year,period,periodicity,lags,tofile);
    local p,pc, lev, kork, aineisto;

    p=acf(res,lags);
    pc=2/sqrt(rows(res));

    begwind;
    lev = 9;
    kork = 6.855;
    makewind(lev,kork/2,0,kork/2,0);
    makewind(lev,kork/2,0,0,1);
    makewind(lev,kork/2,0,0,1);
    /* residual window */
    setwind(1);
    fonts("simplex simgrma");
    ylabel("residuals");
    _pltype={6,6};
    _pcolor=15; /* paksuhko 15 */
    _paxes=1;
    _pnum={1 2};
	if period == 1;
		aineisto = (res)~zeros(rows(res)+period-1,1);
	else;
		aineisto = (crmiss(period-1,1)|res)~zeros(rows(res)+period-1,1);
	endif;
	ts(aineisto,year,periodicity);
    /* autocorrelation window */
    nextwind;
    graphset;
    fonts("simplex simgrma");
    ylabel("\201autocorrelation");
    xlabel("\201lag");
    _pdate="";
    _pbarwid=1-lags/100;
    _pbartyp={6 15};
    _paxes=1;
    _pnum={1 2};
    xtics(0,lags+1,1,0);
    ytics(-1,1,0.1,10);
    bar(seqa(1,1,lags),p);
    /* confidence interval for autocorrelation */
    nextwind;
    _paxes=0;
    _pltype={6 6};
    _pdate="";
    _pcolor=15;
    p=ones(lags+2,2);
    p[.,1]=pc*p[.,1];
    p[.,2]=-pc*p[.,2];
    xy(seqa(0,1,lags+2),p);
    draw;
    if stof(tofile) /= 0;
        graphprt("-C=1 -CS=F -CM=0.1,0.1,0.1,0.1 -CF=" $+ tofile );
    endif;
    endwind;

endp;



/* SKEWNESS
**
** Purpose:   compute skewness of the time series
**
** Format:    y=skewness(x);
**
** Input:     x     NxK matrix, time series
**
** Output:    y     Kx1 matrix, skewnesses
**
** Created:   26.4.94 by AR
** Revised:
**
*/
proc (1)=skewness(x);
     local y,a;

     x=x-meanc(x)';
     y=sumc(x^3)/(rows(x)-1);
     a=(diag(x'x)./(rows(x)-1))^1.5;
     retp(y./a);
endp;


/* KURTOSIS
**
** Purpose:   compute excess kurtosis of the time series
**
** Format:    y=kurtosis(x);
**
** Input:     x     NxK matrix, time series
**
** Output:    y     Kx1 matrix, kurtosis of the series
**
** Created:   26.4.94 by AR
** Revised:
**
*/
proc (1)=kurtosis(x);
     local y,a;

     x=x-meanc(x)';
     y=sumc(x^4)/(rows(x)-1);
     a=(diag(x'x)/(rows(x)-1))^2;
     retp(y./a-3);
endp;



/* JBNORM
**
** Purpose:   compute Jarque Bera normality test statistic, which is
**            chi^2(2) distributed
**
** Format:    y=jbnorm(sk, ek, T, r);
**
** Input:     sk     scalar, skewness
**            ek     scalar, excess kurtosis
**            T      scalar, number of observations
**            r      scalar, number of regressors in regression equation
**
** Output:    y      scalar, the test statistic
**
** Created:   26.4.94 by AR
** Revised:
**
*/
proc (1)=jbnorm(sk, ek, T, r);

     retp((T-r)*(sk^2+ek^2/4)/6);
endp;
