markov.prg

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

FWIG = ZEROS(NGRID,1);  
    
INDEX2 = 1;
do until INDEX2 > NGRID;
Z1 = X0;
Z2 = XGRID;
{PWIG,FWIG} = FWIGGL(Z1,Z2,Z);
INDEX2 = INDEX2 + 1;
endo;

TEMP = MININDC(ABS(FWIG - RN[INDEX]));

if INDEX == Q + 1;
  XTRY = X0|XGRID[TEMP,.];
else;
  XTRY = XTRY|XGRID[TEMP,.];
endif;

N1 = ROWS(XTRY);
N0 = N1 - Q + 1;
X0 = XTRY[N0:N1,.];
INDEX = INDEX + 1;
endo;

retp(XTRY);
endp;

/*******************************************************************************/
/******** ESTIMATION OF H(.) AND T-STATISTIC USING SAMPLE FROM POPULATION ******/
/*******************************************************************************/

@
  HERE WE ESTIMATE THE STATISTIC OF INTEREST USING THE ORIGINAL SAMPLE   
@

proc(2) = HEST(Z,HZERO);

@ HZERO IS THE VALUE OF H UNDER THE NULL HYPOTHESIS @

local MHAT, HHAT, TEMP, VHAT, SHAT, TSTAT, TMP, CY, VY;

MHAT = MEANC(X);  @ MEAN OF VALUES OF X @

@ ***** CALCULATE HHAT, WHICH IS A FUNCTION OF MHAT - WRITE THE FORMULA ****** @
HHAT = MHAT;   @ ***** WRITE FORMULA HERE ****** @
HHAT = HHAT[1];  @ ***** MODIFY INDEX TO DESIRED ESTIMATOR ***** @

@ ****** WRITE THE FORMULA TO FIND VARIANCE OF N^(1/2)*(HHAT - HZERO) ****** @
VY = MEANC(X[.,1]^2) - (MEANC(X[.,1]))^2;
CY = MEANC(Z[.,1].*Z[.,2]);
VHAT = VY + 2*((N - 1)/N)*CY; @ ***** WRITE FORMULA HERE ****** @

@ COMPUTE THE T-STATISTIC @
SHAT = SQRT(VHAT);    @ ***** MODIFY INDEX TO DESIRED ESTIMATOR ***** @
TSTAT = N^(1/2)*(HHAT - HZERO)/SHAT[1];    @ ***** MODIFY INDEX TO DESIRED ESTIMATOR ***** @

retp(HHAT,TSTAT);
endp;

/*******************************************************************************/
/*******************  ESTIMATION of H(.) USING BOOSTRAP SAMPLE *****************/
/*******************************************************************************/

proc(2) = BHEST(ZBOOT);

local BMHAT, BHHAT, BVHAT, BSHAT, XBOOT, BCY, BVY;

XBOOT = ZBOOT[.,1:D];
BMHAT = MEANC(XBOOT);  @ MEAN OF BOOTSTRAP VALUES OF X @

@ ***** CALCULATE BHHAT, WHICH IS A FUNCTION OF MHAT - WRITE THE FORMULA ****** @
BHHAT = BMHAT;   @ ***** WRITE FORMULA HERE ****** @
BHHAT = BHHAT[1];  @ ***** ADJUST INDEX TO DESIRED ESTIMATOR ****** @

@ ****** WRITE THE FORMULA TO FIND VARIANCE OF N^(1/2)*(BHHAT - HHAT) OF EACH BOOTSTRAP 
SAMPLE - THIS FORMULA SHOULD LOOK EXACTLY LIKE THE FORMULA for VHAT, EXCEPT for THE FACT
THAT NOW BOOTSTRAP SAMPLES ARE BEING USED ****** @
BVY = MEANC(X[.,1]^2) - (MEANC(X[.,1]))^2;
BCY = MEANC(Z[.,1].*Z[.,2]);
BVHAT = BVY + 2*((N - 1)/N)*BCY; @ ***** WRITE FORMULA HERE ****** @

@ COMPUTE THE STANDARD DEVIATION @
BSHAT = SQRT(BVHAT[1]);    @ ***** MODIFY INDEX TO DESIRED BOOTSTRAP ESTIMATOR ***** @

retp(BHHAT,BSHAT);
endp;

/*******************************************************************************/
/**************** BOOTSTRAP COMPUTATION OF CRITICAL VALUE **********************/
/*******************************************************************************/

proc(3) = TBOOT(HHAT,Z,ZDAT,NCRIT);

@
  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
@

local XBOOT, INDEX, ZBOOT, BHHAT, BSHAT, IBOOT, BTCUM, BTCRIT1, BTCRIT2,
BSCUM, BHCUM;

IBOOT = 1;
BHCUM = ZEROS(NBOOT,1);  @ STORES BOOTSTRAP ESTIMATES OF H @
BTCUM = ZEROS(NBOOT,1);  @ STORES BOOTSTRAP ESTIMATES OF THE T STATISTIC @
BXCUM = ZEROS(N,D*NBOOT);  @ STORES BOOTSTRAP SAMPLES @
BSCUM = ZEROS(NBOOT,1);  @ STORES BOOTSTRAP ESTIMATES OF THE STDEV @

do until IBOOT > NBOOT;

@ GENERATE A BOOTSTRAP SAMPLE @
XBOOT = BSAMPL(Z,ZDAT);  @ THIS IS N X D @

@ GETTING THE DATA ORGANIZED IN THE APPROPRIATE WAY @
INDEX = 1;
ZBOOT = XBOOT[(Q + 1):N,.];
do until INDEX > Q;
ZBOOT = XBOOT[(Q + 1 - INDEX):(N - INDEX),.]~ZBOOT;
INDEX = INDEX + 1;
endo;

{BHHAT,BSHAT} = BHEST(ZBOOT);
BHCUM[IBOOT] = BHHAT; 
BSCUM[IBOOT] = BSHAT;
BXCUM[.,(1 + (IBOOT - 1)*D):IBOOT*D] = XBOOT;
IBOOT = IBOOT + 1;
endo;

@ CALCULATING THE BOOTSTRAP T-STSTISTIC @
BTCUM = N^(1/2)*(BHCUM - HHAT)./BSCUM;  
BTCUM = SORTC(BTCUM,1);

@ FINDING THE BOOTSTRAP CRITICAL VALUES @
BTCRIT1 = SORTC(ABS(BTCUM),1);  @ TWO - SIDED TEST @
BTCRIT1 = BTCRIT1[NCRIT];       @ TWO - SIDED TEST @
BTCRIT2 = BTCUM[NCRIT];         @ ONE - SIDED TEST @

retp(BTCUM,BTCRIT1,BTCRIT2);
endp;


/***********************************************************************/
/********************* CREATING THE GRID *******************************/
/***********************************************************************/

proc(1) = ONESTEP(A,XGRID);

local INDEX, INDEX1, RK, CK, B;

RK = ROWS(A);
CK = COLS(A);
B = ZEROS(RK*NGRID, CK + 1);
INDEX = 1;
do until INDEX > RK;
INDEX1 = 1;
do until INDEX1 > NGRID;
B[(INDEX1 + (INDEX - 1)*NGRID),.] = A[INDEX,.]~XGRID[INDEX1,(CK + 1)];
INDEX1 = INDEX1 + 1;
endo;
INDEX = INDEX + 1;
endo;

retp(B);
endp;

/***********************************************************************/
/********************* MAIN PROGRAM ************************************/
/***********************************************************************/

TSTART = DATE;

@ ****** MODIFY THE FOLLOWING TO DESIRED VALUES ******* @

NBOOT = 100;   @ ***** NUMBER OF BOOTSTRAP ITERATIONS ***** @
Q = 1;        @ ***** ORDER OF MARKOV PROCESS ***** @
NGRID = 50;   @ ***** NUMBER OF BOOTSTRAP GRID POINTS for EACH X[.,I] - TOTAL NUMBER 
              OF GRID POINTS WILL BE NGRID^D ***** @
LVL = 0.05;   @ ***** NOMINAL LEVEL OF THE TEST ***** @
HX = 3;       @ ***** BANDWIDTH ****** @
HZERO = 0;    @ ***** NULL HYPOTHSIS OF H(.) ****** @

@ ***** CHOOSE LEVEL OF TEST ***** @

TCRIT = {1.96, 1.645};  @ ***** ASYMPTOTIC CRITICAL VALUES for 0.05 LEVEL for 
  2 and 1 SIDED TESTS ****** ({2.575, 2.3264} for 0.01 LEVEL 2 and 1 SIDED) ***** @
NCRIT = CEIL(NBOOT*(1 - LVL));  @ BOOTSTRAP CRITICAL VALUE @


/***************** SETTING UP THE DATA **********************************/

@ 
  load THE DATA and SET UP ARRAYS 
  DATA SHOULD BE AN N X D MATRIX OF OBSERVATIONS OF X, I.E. DATA IS ZDAT
@

load DATASET[50,1] = C:\GAUSS\MARKOV\AR1.TXT; @ ***** SPECIFY SIZE OF MATRIX and PATH TO DATA - DATA SHOULD BE IN A .TXT FILE ***** @
X = DATASET;
N = ROWS(X);  @ NUMBER OF OBSERVATIONS @
D = COLS(X);  @ DIMENSION OF VECTOR X @

@
  STANDARDIZE THE VARIABLES
@
X = (X - MEANC(X)')./STDC(X)';

@ 
  SETUP THE VECTOR Z TO BE USED IN THE MARKOV BOOTSTRAP 
@

INDEX1 = 1;
Z = X[(Q + 1):N,.];
do until INDEX1 > Q;
Z = X[(Q + 1 - INDEX1):(N - INDEX1),.]~Z;
INDEX1 = INDEX1 + 1;
endo;

ZDAT = X;


OUTPUT ON;
"";
"SAMPLE SIZE:  ";; N;
"BOOTSTRAP ITERATIONS:  ";; NBOOT;
"NUMBER OF GRID POINTS:  ";; NGRID;
"ORDER OF MARKOV PROCESS:  ";; Q;
"";
OUTPUT OFF;

/**************** FINDING THE SAMPLE ESTIMATOR OF H(.) ********************/
{HHAT,TSTAT} = HEST(Z,HZERO);

/************** CREATING THE GRID for THE BOOTSTRAP POINTS: XGRID ********/
XMIN = MINC(X);
XMAX = MAXC(X);
XSTEP = ((XMAX) - (XMIN))/(NGRID - 1);
XGRID = ZEROS(NGRID,D);
I = 1;
do until I > D;
XGRID[.,I] = SEQA(XMIN[I],XSTEP[I],NGRID);
I = I + 1;
endo;

if D == 1;
XGRID = XGRID;
else;
A = XGRID[.,1];
II = 1;
do until II > D - 1;
{B} = ONESTEP(A,XGRID);
A = B;
II = II + 1;
endo;
XGRID = B;
endif;

NGRID = NGRID^D;

if NBOOT > 0;

/****** FINDING THE BOOTSTRAP T-STATISTICS and CRITICAL VALUES ************/
{BTCUM,BCRIT1,BCRIT2} = TBOOT(HHAT,Z,ZDAT,NCRIT);

endif;

RTIME = ETHSEC(TSTART,DATE)/6000;  @ RUNNING TIME IN MINUTES @
RTIMEH =ETHSEC(TSTART,DATE)/360000;  @ RUNNING TIME IN HOURS @
OUTPUT ON;
"";
"ASYMPTOTIC CRITICAL VALUES:   ";; TCRIT;
"BOOTSTRAP CRITICAL VALUES:    ";; BCRIT1;; BCRIT2;
"BOOTSTRAP DISTRIBUTION OF T-STATISTIC    ";;  BTCUM;
"RUNNING TIME IN MINUTES    ";; RTIME;
"RUNNING TIME IN HOURS    ";; RTIMEH;
OUTPUT OFF;