/*************************************************************************
**       Gauss Code for "Robust Tests of Extremal Volatility Spillover  **
**                       with an Application to Market Contagion"       **
**                                                                      **
**  Jonathan B. Hill 2007                                               **
**                                                                      **
**  Disclaimer: We make no guarantees that the following code works.    **
**              If an obvious bug exists please send me the details.    **
**              Please do not expect technical support.                 **
*************************************************************************/

new;
cls;
clear all;

library pgraph;

/* It is assumed that two variables in levels are stored in x and y.

   For serial tail dependence, load x, and set y = x.
   
   It is assumed that x and y have the same sample size, and that observations
   represent the same dates. */


    
x=packr(x);             /* default: assumes x, y in levels: transform to log differences */
x=ln(x[.,1]);
x =100.*diff(x,1);

y=packr(y);
y=ln(y[.,1]);
y =100.*diff(y,1);


    /* TEST GLOBALS */
m_pow = .85;                    /* Window of tail fractiles m = [n^m_pow]-[n^m_win]...[n^m_pow]+[n^m_win]  */
m_win = .845;

tail_x = 0;                     /* 0 = two-tailed; 1 = left; 2 = right */
tail_y = 0;

kern_band = .245;               /* kernel var. badnwidth = [m^kern_band] */

h = 1;                          /* Number of horizons  */

tail_coef = 1;                  /* 1 = use tail coef; 2 = use exceed. corr. */

conf_level = .95;                /* for confidence band derivation and plotting */

off_plot = 0;                   /* 1 =  do not plot output */

suppress_out = 1;               /* supresses print otput except final median tail dependence estimates */

/* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */

/* The following ruitine estimates tail spillover over a rolling window of sub-sample sizes,
   starting from some minimum fraction of n, to n;
 
   Set min_sample_fractile < 1 for rolling windows,
   Set min_sample_fractile = 1 for only entire sample. 

   If off_plot = 0, then each window's results will be plotted. 
   If rolling windows are used, it is recommended that suppress_out = 1: only view median estimates and plots*/

z=x~y;
min_sample_fractile = 1;             /* default = 1: entire sample */
increment = .20;                     /* incremental increase in sample fraction used for analysis */

ee=min_sample_fractile;
ee_counter = 1;

do until ee > 1;
l = round(rows(z).*ee);
x=z[1:l,1];y=z[1:l,2];

z_a = cdfni(1-(1-conf_level)/2);

if tail_coef .eq 1;
    {tail_coeff,var_tail_coeff,Wald,CWald,m_tc} = tail_coeff_rolling_m(x,y,h);
    if off_plot .eq 0;
        plot_xy_dep(x,y,tail_coeff,z_a*sqrt(var_tail_coeff./m_tc),title_x,title_y,"Tail Spillover\LRolling Fractile Plots with Robust Intervals");
    endif;
else;
    {exceed_corr,var_exceed_corr,Wald,CWald,m_ec} = exceed_corr_rolling_m(x,y,h);
    if off_plot .eq 0;
        plot_xy_dep(x,y,exceed_corr,z_a*sqrt(var_exceed_corr./m_ec),title_x,title_y,"Tail Spillover\LRolling Fractile with Robust Intervals");
    endif;
endif;

label = "m";
i=1;do until i > h;label = label|"-k";i=i+1;endo;
i=1;do until i > h;label = label|"r";i=i+1;endo;
i=1;do until i > h;label = label|"k";i=i+1;endo;

if tail_coef .eq 1;
    if suppress_out .eq 0;
        $label';
        m_tc~(-z_a*sqrt(var_tail_coeff./m_tc))~tail_coeff~z_a*sqrt(var_tail_coeff./m_tc);"";"";
    endif;
    "                                        MEDIAN TAIL DEPENDENCE VALUES";
    "         h               -k                 r(h)            k                 Wald          p-val";
    seqa(1,1,h)~-median(z_a*sqrt(var_tail_coeff./m_tc))~median(tail_coeff)~median(z_a*sqrt(var_tail_coeff./m_tc))~median(Wald)~cdfchic(median(Wald),seqa(1,1,rows(median(Wald))));
else;
    if suppress_out .eq 0;
        m_ec~(-z_a*sqrt(var_exceed_corr./m_ec))~exceed_corr~z_a*sqrt(var_exceed_corr./m_ec);"";"";
    endif;
    "                                        MEDIAN TAIL DEPENDENCE VALUES";
    "         h               -k                 r(h)            k                 Wald          p-val";
    seqa(1,1,h)~-median(z_a*sqrt(var_exceed_corr./m_ec))~median(exceed_corr)~median(z_a*sqrt(var_exceed_corr./m_ec))~median(Wald)~cdfchic(median(Wald),seqa(1,1,rows(median(Wald))));
endif;

"______________________________________________________________________________________________________";"";"";

"Sample size: " rows(x);


ee=ee+increment;
endo;


/* @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ */


/*------------------------------------------------------------------------
**  tail_coeff_rolling_m.prc
**
**  Purpose: Derive tail dependence coefficient r(x,y,l1,l2,h) 
**           measuring extremal causality from Y(t-h) to X(t)
**
**  Jonathan B. Hill (2006)
-----------------------------------------------------------------------*/

proc (5) = tail_coeff_rolling_m(x,y,h);
    local   i_m,i,j,m,mm,n,n0,tail_coeff;
    local   x0,y0,x0m,y0m,x_b_0,y_b_0,x0_lag,y0_lag;
    local   z_b_0,z_b_mat,m_min,m_max,var;
    local   z_b_0_matrix,cov_h,m_cov_h,W,CW;

    if tail_x .eq 0;
        x0 = abs(x);
    endif;
    if tail_x .eq 1;
        x0 = -x.*(x .< 0);
    endif;
    if tail_x .eq 2;
        x0 = x.*(x .> 0);
    endif;
    if tail_y .eq 0;
        y0 = abs(y);
    endif;
    if tail_y .eq 1;
       y0 = -y.*(y .< 0);
    endif;
    if tail_y .eq 2;
        y0 = y.*(y .> 0);
    endif;

    x0_lag = packr(lagn(x0,seqa(0,1,h+1)));
    y0_lag = packr(lagn(y0,seqa(0,1,h+1)));

    n = rows(x0_lag);

    n0 = minc(sumc(x0 .> 0)|sumc(y0 .> 0));

    m_min = round(n0^m_pow)-round(n0^m_win);
    m_max = round(n0^m_pow)+round(n0^m_win);

    tail_coeff=zeros(m_max-m_min+1,h);
    var=zeros(m_max-m_min+1,h);

    mm=seqa(m_min,1,m_max-m_min+1);

    z_b_0_matrix = zeros(rows(x0_lag),h);

    W = zeros(m_max-m_min+1,h);
    CW = zeros(m_max-m_min+1,h);

    i_m = 1;
    do until i_m > rows(mm);

        m = mm[i_m,1];

        z_b_mat = zeros(rows(x0_lag),h);

        i=1;
        do until i >h;

            x0m=sortc(x0_lag[.,1],1);
            y0m=sortc(y0_lag[.,i+1],1);

            x_b_0 = (ln(x0_lag[.,1]) .> ln(x0m[n-m+1,1]));
            y_b_0 = (ln(y0_lag[.,i+1]) .> ln(y0m[n-m+1,1]));

            z_b_0 = packr(x_b_0~y_b_0);
            x_b_0 = z_b_0[.,1];
            y_b_0 = z_b_0[.,2];

            z_b_0 = (x_b_0-(m/n)).*(y_b_0-(m/n));
            z_b_0_matrix[.,i] = z_b_0;

            tail_coeff[m-m_min+1,i] = sumc(z_b_0)/m;

            {var[m-m_min+1,i]} = kern_var_general(z_b_0,m);

            i=i+1;

        endo;

        {cov_h,m_cov_h} = kern_cov_matrix_general(z_b_0_matrix,m);

        i=1;
        do until i > h;
            W[m-m_min+1,i] = (sqrt(m_cov_h[1:i,1]').*tail_coeff[m-m_min+1,1:i])*inv(cov_h[1:i,1:i])*(sqrt(m_cov_h[1:i,1]).*tail_coeff[m-m_min+1,1:i]');
            CW[m-m_min+1,i] = ((W[m-m_min+1,i]-i)^2)/(2*i);
            i=i+1;
        endo;

        i_m=i_m+1;
    endo;

    retp(tail_coeff,var,W,CW,mm);
endp;
/*____________________________________________________________________________________*/


/*------------------------------------------------------------------------
**  exceed_corr_rolling_m.prc
**
**
**  Jonathan B. Hill (2006)
-----------------------------------------------------------------------*/

proc (5) = exceed_corr_rolling_m(x,y,h);
    local   i_m,i,j,m,mm,n,n0,exceed_corr;
    local   x0,y0,x0m,y0m,x_b_0,y_b_0,x0_lag,y0_lag;
    local   z_b_0,m_min,m_max,s_x_b_0,s_y_b_0;
    local   a_hat_x,a_hat_y,cov,cov_a_inv,var;
    local   z_b_0_matrix,W,CW,cov_h,m_cov_h;

    if tail_x .eq 0;
        x0 = abs(x);
    endif;
    if tail_x .eq 1;
        x0 = -x.*(x .< 0);
    endif;
    if tail_x .eq 2;
        x0 = x.*(x .> 0);
    endif;
    if tail_y .eq 0;
        y0 = abs(y);
    endif;
    if tail_y .eq 1;
       y0 = -y.*(y .< 0);
    endif;
    if tail_y .eq 2;
        y0 = y.*(y .> 0);
    endif;

    x0_lag = packr(lagn(x0,seqa(0,1,h+1)));
    y0_lag = packr(lagn(y0,seqa(0,1,h+1)));

    n = rows(x0_lag);

    n0 = minc(sumc(x0 .> 0)|sumc(y0 .> 0));

    m_min = round(n0^m_pow)-round(n0^m_win);
    m_max = round(n0^m_pow)+round(n0^m_win);

    exceed_corr=zeros(m_max-m_min+1,h);
    var = zeros(m_max-m_min+1,h);

    mm=seqa(m_min,1,m_max-m_min+1);

    z_b_0_matrix = zeros(rows(x0_lag),h);

    W = zeros(m_max-m_min+1,h);
    CW = zeros(m_max-m_min+1,h);

    i_m = 1;
    do until i_m > rows(mm);

        m = mm[i_m,1];

        i=1;
        do until i >h;

            x0m=sortc(x0_lag[.,1],1);
            y0m=sortc(y0_lag[.,i+1],1);

            x_b_0 = (ln(x0_lag[.,1]) - ln(x0m[n-m+1,1])).*(ln(x0_lag[.,1]) .> ln(x0m[n-m+1,1]));
            y_b_0 = (ln(y0_lag[.,i+1]) - ln(y0m[n-m+1,1])).*(ln(y0_lag[.,i+1]) .> ln(y0m[n-m+1,1]));

            z_b_0 = packr(x_b_0~y_b_0);
            x_b_0 = z_b_0[.,1];
            y_b_0 = z_b_0[.,2];

            {a_hat_x,a_hat_y,cov,cov_a_inv} = a_hat_kern_var_xy(x,y,m);

            z_b_0 = (x_b_0-(m/n)*(1./a_hat_x)).*(y_b_0-(m/n)*(1./a_hat_y))./(2*(1./a_hat_x)*(1./a_hat_y));
            z_b_0_matrix[.,i] = z_b_0;
            exceed_corr[m-m_min+1,i] = sumc(z_b_0)/m;

            {var[m-m_min+1,i]} = kern_var_general(z_b_0,m);

            i=i+1;

        endo;

        {cov_h,m_cov_h} = kern_cov_matrix_general(z_b_0_matrix,m);

        i=1;
        do until i > h;
            W[m-m_min+1,i] = (sqrt(m_cov_h[1:i,1]').*exceed_corr[m-m_min+1,1:i])*inv(cov_h[1:i,1:i])*(sqrt(m_cov_h[1:i,1]).*exceed_corr[m-m_min+1,1:i]');
            CW[m-m_min+1,i] = ((W[m-m_min+1,i]-i)^2)/(2*i);
            i=i+1;
        endo;

        i_m=i_m+1;
    endo;

    retp(exceed_corr,var,W,CW,mm);
endp;
/*____________________________________________________________________________________*/


/*-------------------------------------------------------------------------
**  kern_var_general.prc
**
**  Purpose:        Compute kernel variance estimator for some array {z(n,t)}
**                  
**  Usage:          {var} = kern_var_general(z,m)
**
**  Input:          z       n-vector
**                  m       sample fractile
**
**  Output          var     kernel variance estimator
**
**  Comments:       Requires global "kern_band"
**
**  Jonathan B. Hill (2006)
-----------------------------------------------------------------------*/

proc (1) = kern_var_general(z,m);
    local   i,j,k,w,n,bn,z_lag,var;

    /*z=abs(z);*/

    z=packr(miss(z,0));

    n = rows(z);

    i=1;
    do until i > rows(m);
        j=1;
        do until j > cols(m);

            if m[i,j] .eq 0;
                goto skip_kern_var_general;      /* reset m = 1 if ranked m is simply "0" */
            endif;

            bn = round(m[i,j].^kern_band);
            var = sumc(z[.,1].^2)./m[i,j];       /* asymptotic variance */

            k=1;
            do until k > bn;

                z_lag = packr(lagn(z,(0|k)));
                w = 1-k/bn;                     /* Bartlett kernel */
                var = var + 2*w*sumc(z_lag[.,1].*z_lag[.,2])./m[i,j];
                k=k+1;
            endo;

            skip_kern_var_general:
            j=j+1;
        endo;
        i=i+1;
    endo;

    retp(var);

endp;

/*____________________________________________________________________________________*/



/*-------------------------------------------------------------------------
**  kern_cov_matrix_general.prc
**
**  Purpose:  Derive kernal asymp. cov. mat. estimator for h-vector of tail dependence estimates
**                  
**  Usage:     {cov} = a_hat_kern_var_z_vec(zz,m);
**
**  Input:      zz      nxh time series
**              m       sample tail fractile
**
**  Output      cov     un-scaled: cov and not cov/m
**              mz              
**
**  Jonathan B. Hill (2006)
-----------------------------------------------------------------------*/

proc (2) = kern_cov_matrix_general(zz,m);
    local   nz,nzz,mx,mz,mzz,i,j,h,w;
    local   cov,n,z_ord,n_pos,a_inv_z,bn,ln_zz;
    local   zz_lag_mj,zzs,zz_lag_pj;

    n = rows(zz);h = cols(zz);nz = sumc(zz .> 0);
    mz = round(nz^(ln(m)/ln(n)));
    mzz = minc(mz);
    nzz = minc(nz);
    cov = miss(zeros(h,h),0);
    bn = round(mzz.^kern_band);

    cov = zz'zz/mzz;

    j=1;
    do until j > bn;
        zz_lag_mj = zeros(rows(zz),h);
        zz_lag_pj = zeros(rows(zz),h);
        zzs=zz;

        i=1;
        do until i > h;
            zz_lag_mj[.,i] = lagn(zzs[.,i],j);
            zz_lag_pj[.,i] = lagn(zzs[.,i],-j);
            i=i+1;
        endo;

        zz_lag_pj=packr(zzs~zz_lag_mj~zz_lag_pj);
        zzs = zz_lag_pj[.,1:cols(zzs)];
        zz_lag_mj = zz_lag_pj[.,cols(zzs)+1:cols(zzs)+cols(zz_lag_mj)];
        zz_lag_pj = zz_lag_pj[.,cols(zzs)+cols(zz_lag_mj)+1:cols(zz_lag_pj)];

        w = 1-j/bn;

        cov = cov + w*(zzs'zz_lag_mj+zz_lag_pj'zzs)./mzz;

        j=j+1;
    endo;

    retp(cov,mz);

endp;

/*____________________________________________________________________________________*/




/*-------------------------------------------------------------------------
**  a_hat_kern_var_xy.prc
**
**  Purpose:  Derive Hill's (1975) tail index estimator for two processes {x(t),y(t)}
**            with support on [0,infin), and a kernel estimator of the covariance matrix, 
**            cf. J.B. Hill (2006)
**                  
**  Usage:     {a,cov} = a_hat_kern_var_xy(x,m);
**
**  Input:      x       nx1 time series
**              y       nx1 time series
**              m       order index used for a,c derivation: assumed that m is a scalar
**
**  Output      a_hat_x         x index estimator
**              a_hat_y         y index estimators
**              cov             kernel variance estimator
**              cov_a_inv       kernel variance estimator of 1/a
**
**  Jonathan B. Hill (2006)
-----------------------------------------------------------------------*/

proc (4) = a_hat_kern_var_xy(x,y,m);
    local   nx,ny,nxy,mx,my,mxy,i,j,k,w,a_hat_x,a_hat_y,z;
    local   cov,n,x_ord,y_ord,n_pos;
    local   a_inv_x,a_inv_y,bn,ln_xx,ln_yy;
    local   ln_xx_lag,ln_yy_lag,cov_a_inv,tail_flag,var;

    n = rows(x);

    nx=sumc(x .> 0);
    ny=sumc(y .> 0);

    cov_a_inv = miss(zeros(2,2),0);
    cov = miss(zeros(2,2),0);

    if nx .> 0;
        tail_flag = 1;
        if ny .> 0;
            tail_flag = 12;
            mx = round(nx^(ln(m)/ln(n)));
            my = round(ny^(ln(m)/ln(n)));
            mxy = minc(mx|my);              /* DOES THIS MATCH OTHER m12 settings? */
            nxy = minc(nx|ny);
        endif;
    else;
        if ny .> 0;
            tail_flag = 2;
        else;
            tail_flag = 0;
        endif;
    endif;

    if tail_flag .eq 0;
        "Tail Index Estimation ERROR: X(t) = 0 for all t";
        retp(0,0,0,0);
    elseif tail_flag .eq 1;
        {a_hat_x,var} = a_hat_kern_var(x,m);
        a_hat_y = 0;
        cov[1,1] = var;
        cov_a_inv[1,1] = cov[1,1].*(a_hat_x^(-2));
    elseif tail_flag .eq 2;
        {a_hat_y,var} = a_hat_kern_var(y,m);
        a_hat_x = 0;
        cov[2,2] = var;
        cov_a_inv[2,2] = cov[2,2].*(a_hat_y^(-2));
    else;

        x_ord = sortc(x,1);
        x_ord = x_ord[n-nxy+1:n,1];

        y_ord = sortc(y,1);
        y_ord = y_ord[n-nxy+1:n,1];

        a_hat_x = (meanc(ln(abs(x_ord[nxy-mxy+1:nxy,1])))-ln(abs(x_ord[nxy-mxy,1])))^(-1);
        a_hat_y = (meanc(ln(abs(y_ord[nxy-mxy+1:nxy,1])))-ln(abs(y_ord[nxy-mxy,1])))^(-1);

            /* compute newey west variance estimator for a_hat^(-1) for het., dep. data */

        a_inv_x = 1./a_hat_x;
        a_inv_y = 1./a_hat_y;

        bn = round(minc(mx|my).^kern_band);

        ln_xx = ((ln(x_ord[.,1])-ln(x_ord[nxy-mxy,1])).*(ln(x_ord[.,1])-ln(x_ord[nxy-mxy,1]) .> 0) - (mxy/nxy)*a_inv_x);
        ln_yy = ((ln(y_ord[.,1])-ln(y_ord[nxy-mxy,1])).*(ln(y_ord[.,1])-ln(y_ord[nxy-mxy,1]) .> 0) - (mxy/nxy)*a_inv_y);
        cov_a_inv[1,1] = sumc(ln_xx[.,1].^2)./mxy;
        cov_a_inv[1,2] = sumc(ln_xx[.,1].*ln_yy[.,1])./mxy;
        cov_a_inv[2,1] = sumc(ln_xx[.,1].*ln_yy[.,1])./mxy;
        cov_a_inv[2,2] = sumc(ln_yy[.,1].^2)./mxy;

        k=1;
        do until k > bn;

            ln_xx_lag = packr(lagn(ln_xx,(0|k)));
            ln_yy_lag = packr(lagn(ln_yy,(0|k)));
            w = 1-k/bn;                               /* Bartlett kernel */
            cov_a_inv[1,1] = cov_a_inv[1,1]  + 2*w*sumc(ln_xx_lag[.,1].*ln_xx_lag[.,2])./mxy;
            cov_a_inv[1,2] = cov_a_inv[1,2]  - 2*w*sumc(ln_xx_lag[.,1].*ln_yy_lag[.,2])./mxy;
            cov_a_inv[2,1] = cov_a_inv[2,1]  - 2*w*sumc(ln_yy_lag[.,1].*ln_xx_lag[.,2])./mxy;     
            cov_a_inv[2,2] = cov_a_inv[2,2]  + 2*w*sumc(ln_yy_lag[.,1].*ln_yy_lag[.,2])./mxy;

            k=k+1;
        endo;

        cov = abs((cov_a_inv).*((a_hat_x^2|a_hat_y^2)*(a_hat_x^2|a_hat_y^2)'));

            /* scale as asymptotic covariances */
        cov[1,1]=cov[1,1]/mx;cov[2,2]=cov[2,2]/my;cov[1,2]=cov[1,2]/sqrt(mx*my);cov[2,1]=cov[1,2];
        cov_a_inv[1,1]=cov_a_inv[1,1]/mx;cov_a_inv[2,2]=cov_a_inv[2,2]/my;cov_a_inv[1,2]=cov_a_inv[1,2]/sqrt(mx*my);cov_a_inv[2,1]=cov[1,2]/sqrt(mx*my);

        cov[1,2] = (cov[1,2]+cov[2,1])/2;cov[2,1] = cov[1,2];
        cov_a_inv[1,2] = (cov_a_inv[1,2]+cov_a_inv[2,1])/2;cov_a_inv[2,1] = cov_a_inv[1,2];

    endif;

    retp(a_hat_x,a_hat_y,cov,cov_a_inv);

endp;
/*____________________________________________________________________________________*/


/*-------------------------------------------------------------------------
**  a_hat_kern_var.prc
**
**  Purpose:  Derive Hill's (1975) tail index estimator a, and Newey-West variance estimator
**                  
**  Usage:     {a_hat,var} = a_hat_kern_var(x,m);
**
**  Input:      x       nx1 time series
**              m       order index used for a,c derivation: m = o(n^d) for some 0 < d < 1
**
**  Output      a_hat       tail index estimator
**              var         kernel variance estimator
**
**  Jonathan B. Hill (2006)
-----------------------------------------------------------------------*/

proc (2) = a_hat_kern_var(x1,m);
    local   i,j,k,w,a_hat,z,x_m_a,n,x_ord,n_pos;
    local   a_inv,bn,ln_xx,ln_xx_lag,var,var_a;

    x1=abs(x1);             /* a_hat for all x: |x| */
    x1=packr(miss(x1,0));   /* effectively removes "0's", and does not change true "alpha" */

    n = rows(x1);

    x_ord = sortc(x1,1);

    a_hat = zeros(rows(m),cols(m));
    var_a = zeros(rows(m),cols(m));

    i=1;
    do until i > rows(m);
        j=1;
        do until j > cols(m);

            if m[i,j] .eq 0;
                goto skip_a;         /* reset m = 1 if ranked m is simply "0" */
            endif;

            a_hat[i,j] = (meanc(ln(abs(x_ord[n-m[i,j]+1:n,1])))-ln(abs(x_ord[n-m[i,j],1])))^(-1);

            a_inv = 1./a_hat[i,j];
            bn = round(m[i,j].^kern_band);

            ln_xx = (ln(x1[.,1])-ln(abs(x_ord[n-m[i,j],1])));
            x_m_a = x_ord[n-m[i,j],1].*a_hat[i,j];              /* X(m+1)^a_hat */
            ln_xx = ln_xx.*(ln_xx .> 0) - (m[i,j]/n)*a_inv;     /* used to estimate asymptotic variance */

            var = sumc(ln_xx[.,1].^2)./m[i,j];       /* asymptotic variance */

            k=1;
            do until k > bn;
                ln_xx_lag = packr(lagn(ln_xx,(0|k)));
                w = 1-k/bn;             /* Bartlett kernel */
                var = var + 2*w*sumc(ln_xx_lag[.,1].*ln_xx_lag[.,2])./m[i,j];
                k=k+1;
            endo;

            var = abs((var).*a_hat[i,j]^4);

            var_a[i,j] = var./m[i,j];      /* scaled for interval estimation*/

            skip_a:
            j=j+1;
        endo;
        i=i+1;
    endo;

    retp(a_hat,var_a);

endp;

/*____________________________________________________________________________________*/


/*--------------------------------------------------------------------------
**  diff.prc
----------------------------------------------------------------------------*/

proc diff(W,d);
    local   i,W1,W2;

    if d == 0;
        retp(W);
    endif;
    if round(d) .ne d;
        "Error: difference lag is not an integer";
        retp(error(0));
    else;
        if d .< 0;
            "Error: difference lag < 0";
             retp(error(0));
        endif;
    endif;    
    W1 = zeros(rows(W)-d,cols(W));
    W2 = zeros(rows(W),d+1);      
    i=1;
    do until i > cols(W);
        W2 = packr(W[.,i]~lagn(W[.,i],d));
        W1[.,i] = W2[.,1]-W2[.,2];
        i=i+1;
    endo;
    
    retp(W1);
endp;

/*____________________________________________________________________________________*/


proc (0) = plot_xy_dep(x,y,rho,k,title_x,title_y,title_dep);
    local   max_y;

    begwind;
    window (3,1,1);
    setwind(1);
  
    graphset;
    _pmcolor =0;
    _ptitlht = .25;_pnumht = .25;_paxht = .25;
    title (title_x);
    _plegctl = 0;_plegstr = "";
    xlabel("day");ylabel(" ");
    xy(0,x);
    nextwind;

    graphset;
    _pmcolor =0;
    _ptitlht = .25;_pnumht = .25;_paxht = .25;
    title (title_y);
    _plegctl = 0;_plegstr = "";
    xlabel("day");ylabel(" ");
    xy(0,y);
    nextwind;

    graphset;
    _pmcolor =0;
    _ptitlht = .25;_pnumht = .25;_paxht = .25;
    title (title_dep);
    _plegctl = 0;_plegstr = "";
    xlabel("m");ylabel(" ");
    max_y = round(100.*maxc(rho|k))./100;
    xy(0,-k~rho~k);
    graphset;
    endwind; 
endp;

/*____________________________________________________________________________________*/