markov.prg

一个重要的关于马尔科夫转移概率的计量经济应用程序

/*************************************************************/ 
/*************************************************************/
/******************* MARKOV **********************************/
/*************************************************************/
/***** PROGRAM TO CARRY OUT BOOTSTRAP ESTIMATION *************/
/***** OF CONFIDENCE INTERVALS OF AUTOREGRESSIVE *************/
/***** MARKOV PROCESSES **************************************/
/*************************************************************/
/*************************************************************/
/****** 27 SEPTEMBER 2004 ************************************/
/*************************************************************/

/*************************************************************/
/** REFERENCE: BOOTSTRAP METHODS for MARKOV PROCESSES (2003) */
/** ECONOMETRICA, VOL. 71, NO.4 ******************************/


new;
FORMAT 9, 4;
clear N, D, X, Z, HX, Q, DELX, XGRID, NGRID, NBOOT, LVL, BXCUM, A, B;
OUTPUT FILE = C:\GAUSS\MARKOV\MARKOV1.OUT ON;

@ *********************** NOTATION *********************************  @
@ N = SAMPLE SIZE
  D = NUMBER OF COLUMS OF X
  X = ORIGINAL DATA
  Z = DATA VECTOR OF X and Y
  HX = BANDWIDTH
  Q = ORDER OF MARKOV PROCESS
  DELX = WIDTH OF GRIDCELLS for BOOTSTRAP VALUES OF X - D X 1 VECTOR
  NG = NUMBER OF GRID POINTS IN EACH DIMENSION OF XGRID
  XGRID = GRID OF X POINTS for BOOTSTRAP VALUES OF X - NG^D MATRIX
  NGRID = NUMBER OF ELEMENTS IN XGRID
  NBOOT = NUMBER OF BOOTSTRAP SAMPLES
  LVL = NOMINAL LEVEL OF A TEST
@

/******************************************************************************/
/*********************SECOND ORDER KERNEL WITH SUPPORT [-1, 1] ****************/
/******************************************************************************/

@
  GIVES THE SECOND ORDER KERNEL for ARRAY V WITH SUPPORT [-1, 1] 
@

proc(1) = KERN2(V,K);

@
  if K < 1.5, COMPUTE KERNEL. OTHERWISE, COMPUTE INTEGRAL OF KERNEL
@ 

local TERMTEMP, TERM, TERM1, TERM2, TERM3, TERM4, VTERM4, INDEX, INDEX2, INDEX3, INDEX4, INDEX5,
TMP, TMP2, TMP3, TMP4, TMP5, VTM4, NV, CV, VTM4MEAN, I;

NV = ROWS(V);
CV = COLS(V);
TERM = ZEROS(NV,1);
I = 1;
if K < 1.5;
  TERM = ZEROS(NV,1);
  do until I > NV;
  TERM[I] = (15/16).*((1 - V[I,.]*V[I,.]')^2);
  TERM[I] = TERM[I].*(V[I,.]*V[I,.]' .<= 1);
  I = I + 1;
  endo;
else;
  if CV == 1;
    TERM = ZEROS(NV,1);
    TERM = (15/16)*((1/5)*V.^5 + V - (2/3)*V.^3) + 0.5;
    TERM = TERM.*(ABS(V) .<= 1) + (V .> 1);
  else;
    TMP = V[.,1];
    INDEX = 2;
    do until INDEX > CV;
    TMP = TMP.*V[.,INDEX];
    INDEX = INDEX + 1;
    endo;
    TERM1 = TMP;
    
    TERM2 = V.^4;
    TERM3 = V.^2;
    TERM2 = TERM2.*TERM1;
    TERM3 = TERM3.*TERM1;
    TMP2 = TERM2[.,1];
    TMP3 = TERM3[.,1];
    INDEX2 = 2;
    do until INDEX2 > CV;
    TMP2 = TMP2 + TERM2[.,INDEX2];
    TMP3 = TMP3 + TERM3[.,INDEX2];
    INDEX2 = INDEX2 + 1;
    endo;
    TERM2 = TMP2;
    TERM3 = TMP3; 

    VTERM4 = ZEROS(NV,CV - 1);
    INDEX3 = 1;
    do until INDEX3 > CV - 1;
    VTM4 = ZEROS(NV,CV - 1);
    INDEX4 = INDEX3 + 1;
    do until INDEX4 > CV;
    TMP4 = (V[.,INDEX3].^2).*(V[.,INDEX4].^2);
    TMP4 = TMP4.*TERM1;
    VTM4[.,INDEX3] = TMP4;
    INDEX4 = INDEX4 + 1;
    endo;
    
    VTM4 = SUMC(VTM4');
    VTERM4[.,INDEX3] = VTM4;
    INDEX3 = INDEX3 + 1;
    endo;
    TERM4 = SUMC(VTERM4');
    
    TERM = ZEROS(NV,1);
    TERM = (15/16)*(TERM1 + (1/5)*TERM2 - (2/3)*TERM3 + (2/9)*TERM4) + 0.5;
    do until I > NV;
    TERM[I] = TERM[I].*(V[I,.]*V[I,.]' .<= 1);
    I = I + 1;
    endo;
  endif;
   
endif;

retp(TERM);
endp;

/******************************************************************************/
/********* COMPUTE CONDITIONAL CDF AND ALL NECESSARY DENSITIES ****************/
/******************************************************************************/

proc(2) = FWIGGL(Z1,Z2,Z);

@
  THIS PROCEDURE COMPUTES DENSITY ESTIMATORS GWIGGL, PWIGGL and FWIGGL
  (CONDITIONAL CDF). GWIGGL IS COMPUTED AT POINT Z1. PWIGGL IS COMPUTED
  AT POINTS X[Q+1],...,X[2] and X[Q],...,X[1], WHERE X[.] IS A D-VECTOR. 
  GWIGGL, PWIGGL and FWIGGL ARE NGRID X 1.
@

@
  Z = (N - Q) X D*(Q + 1): MATRIX CONTAINING OBSERVATIONS OF THE STATE
  VARIABLE X and ITS Q LAGGED VALUES, Y. Z IS THE DATA. IMPORTANT: THE FIRST ROWS 
  OF Z SHOULD HAVE THE LOWEST TIME INDEXES, AS WELL AS THE FIRST COLUMNS. 
  THEREFORE Z = X[1:N - Q,.],...,X[Q + 1:N,.] 
  

  Z1 = Q X D VECTOR CONSISTING OF CURRENT and LAGGED VALUES OF THE 
  STATE VARIABLE. THE FIRST ROW OF Z1 SHOULD HAVE THE LOWEST TIME INDEX.

  Z2 = NGRID X D VECTOR OF POSSIBLE VALUES OF STATE VARIABLE IN THE
  NEXT PERIOD. THE TRANSITION PROBABILITY IS OF Z2 CONDTIONAL ON Z1.

@

local INDEX, NV, J, TEMP, DIF, INDEX1, DZTEMP, DYPTEMP, FWIG, PPWIG, GWIG, PWIG,
INDEX2, INDEX3, INDEX4, INDEX5, DYP, DY, DZZ, DZ;

INDEX = 1;
J = D*(INDEX - 1);
do until INDEX > Q;
if INDEX == 1;
  DZ = KERN2((Z[.,(1 + J):(J + D)] - Z1[INDEX,.])/HX,1);
  DY = DZ;
  DYP = ONES(N - Q,1);
else;
  DZ = DZ.*KERN2((Z[.,(1 + J):(J + D)] - Z1[INDEX,.])/HX,1);
  DY = DZ;
  DYP = DYP.*KERN2((Z[.,(1 + J):(J + D)] - Z1[INDEX,.])/HX,1);
endif;
INDEX = INDEX + 1;
endo;

DIF = ZEROS((N - Q)*NGRID,D);
INDEX1 = 1;
do until INDEX1 > (N - Q);
DIF[(1 + (INDEX1 - 1)*NGRID):(INDEX1*NGRID),.] = Z[INDEX1,(D*Q + 1):(D*(Q + 1))] - Z2;
INDEX1 = INDEX1 + 1;
endo;

DZTEMP = ZEROS((N - Q)*NGRID,1);
DYPTEMP = ZEROS((N - Q)*NGRID,1);
INDEX2 = 1;
do until INDEX2 > (N - Q);
INDEX3 = 1;
do until INDEX3 > NGRID;
DZTEMP[(INDEX3 + (INDEX2 - 1)*NGRID),.] = DZ[INDEX2,.];
DYPTEMP[(INDEX3 + (INDEX2 - 1)*NGRID),.] = DYP[INDEX2,.];
INDEX3 = INDEX3 + 1;
endo;
INDEX2 = INDEX2 + 1;
endo;

DZZ = DZTEMP.*KERN2(DIF/HX,2); @ CDF OF Z[.,1:Q] @
DZ = DZTEMP.*KERN2(DIF/HX,1);  @ DENSITY OF Z[.,1],...,Z[.,Q + 1] @
DYP = DYPTEMP.*KERN2(DIF/HX,1); @ DENSITY OF Z[.,2],...,Z[.,Q + 1] @

GWIG = ZEROS(NGRID,N - Q);
FWIG = ZEROS(NGRID,N - Q);
PPWIG = ZEROS(NGRID,N - Q);
INDEX4 = 1;
do until INDEX4 > NGRID; 
INDEX5 = 1;
do until INDEX5 > (N - Q);
GWIG[INDEX4,INDEX5] = DZ[(INDEX4 + (INDEX5 - 1)*NGRID),.];
FWIG[INDEX4,INDEX5] = DZZ[(INDEX4 + (INDEX5 - 1)*NGRID),.];
PPWIG[INDEX4,INDEX5] = DYP[(INDEX4 + (INDEX5 - 1)*NGRID),.];
INDEX5 = INDEX5 + 1;
endo;
INDEX4 = INDEX4 + 1;
endo;

GWIG = MEANC(GWIG');
FWIG = MEANC(FWIG');
PPWIG = MEANC(PPWIG');
PWIG = MEANC(DY);

GWIG = GWIG/(HX^(Q + 1));  @ DENSITY OF Z[.,1],...,Z[.,Q + 1] @
FWIG = FWIG/(HX^Q);  @ CDF AT Z1[Q + 1], DENSITY OF Z[.,1:Q] @
PWIG = PWIG/(HX^Q); @ DENSITY OF Z[.,1],...,Z[.,Q] @
TEMP = PWIG.*(ABS(PWIG) .> 1e-10) + 1e-10*(ABS(PWIG) .<= 1e-10);
PPWIG = PPWIG/(HX^Q);  @ DENSITY OF Z[.,2],...,Z[.,Q + 1] @
FWIG = FWIG./TEMP;

retp(PPWIG,FWIG); @ PPWIG IS THE MARGINAL DENSITY and FWIG IS THE TRANSITION CDF @
endp;

/***********************************************************************/
/************** GENERATE INITIAL BOOTSTRAP OBSERVATION *****************/
/***********************************************************************/

@
  THIS PROCEDURE GENERATES THE INITIAL STATE OF THE BOOTSTRAP MARKOV
  PROCESS BY SAMPLING Q CONSECUTIVE OBSERVATIONS FROM THE ORIGINAL DATA

  ZDAT = N X D MATRIX CONTAINING THE DATA
  Z[1,.] HAS THE LOWEST TIME INDEX
@

proc(1) = XSTART(ZDAT);

local N1, N2, ZVEC;

N1 = CEIL((N - Q + 1)*RNDU(1,1));
N2 = N1 + Q - 1;
ZVEC = ZDAT[N1:N2,.];

retp(ZVEC);
endp;

/***********************************************************************/
/************** GENERATE BOOTSTRAP SAMPLE ******************************/
/***********************************************************************/

proc(1) = BSAMPL(Z,ZDAT);

@
  THIS PROCEDURE GENERATES ALL MARKOV BOOTSTRAP OBSERVATIONS EXCEPT THE 
  FIRST. THE FIRST IS GENERATED IN THE PROCEDURE XSTART

  THE OUTPUT IS A COLUMN VECTOR OF N OBSERVATIONS OF A BOOSTRAP SCALAR X
  
  AS BEFORE, Z = (N - Q) X (D*Q + D): MATRIX CONTAINING OBSERVATIONS OF 
  THE STATE VARIABLE X and ITS Q LAGGED VALUES, Y. Z IS THE DATA. IMPORTANT: 
  THE FIRST ROWS OF Z SHOULD HAVE THE LOWEST TIME INDEXES, AS WELL AS THE 
  FIRST COLUMNS. THEREFORE Z = X[1:N - Q,.],...,X[Q + 1:N,.] 
   
  ZDAT = N X D MATRIX OF ORIGINAL DATA
  
  OUTPUT IS A SERIES OF N BOOTSTRAP OBSERVATIONS
@

local INDEX, Z1, INDEX2, PWIG, FWIG, RN, TEMP, XTRY, N0, N1, Z2, XTMP, X0;

X0 = XSTART(ZDAT);
RN = RNDU(N,1);

INDEX = Q + 1;
do until INDEX > N;