* * CYCLES.RAT * *************************************************************** * * * * This program calculates business cycle statistics using * * 1. HP-filtered data, * 2. The band-pass filters discussed in Baxter and King, * (this filter includes all frequencies with a period * of less than X years) * 3. The raw data, * 4. First-differenced data. * * It calculates the standard errors of these estimates * * using * * 1. The VARHAC procedure from Den Haan and Levin (1995), * 2. The optimal bandwidth procedure from Newey & West (1994). * * * The user has to adjust the "READING THE DATA" section, * and the "PARAMETER" section below. * * ******************************************************************** ******************************************************************** * PARAMETERS ******************************************************************** ******************************************************************** calendar 1900 1 4 allocate 1999:4 compute [integer] nfirst = 1954:1 compute [integer] nlast = 1994:2 compute [integer] ntot = nlast - nfirst + 1 compute [integer] nvar = 5 compute [integer] start = 1954:1 compute [integer] end = 1994:2 compute [integer] ihp = 1 compute [integer] ibk = 1 compute [integer] iraw = 0 compute [integer] idif = 0 compute [integer] ivarhac = 1 compute [integer] inw = 1 compute [integer] bkper = 32 compute [integer] bklost = 12 * * nfirst: date of the first observation in data file * nlast: date of last observation in data file * ntot: total number of observations in the data file * nvar: total number of variables in the data file * start: first observation to be used ( start.ge.nfirst) * end: last observation to be used ( end.le.nlast ) * ihp: if 1, then statistics for hp-filter will be calculated * ibf: if 1, then statistics for bandpass filters will be calculated * iraw: if 1, then statistics for raw series are calculated * idif: if 1, then statistics for first-differences will be calculated * ivarhac: if 1, then VARHAC standard errors will be calculated * inw: if 1, then NW standard errors will be calculated * bkper: maximum period of the included cycles. Thus, if bkper = 32, * then that part of the series with cycles less than or equal * to 32 periods, i.e. 8 years, is included. * bklost: the BK-filter is a two-sided moving average. The higher * the order the more accurate the filter but the more observations * are lost. If start = 1 and end = ntot, then the statistics * are calculated over the period from start+bklost to end-bklost. * ******************************************************************** ******************************************************************** * READING THE DATA ******************************************************************** ******************************************************************** declare rectangle rawdat(ntot,nvar) open data cycles.dat read rawdat open copy cycles.out * * the user only has to adjust the name of the input and output file, * this program assumes that data are by observation in the file * cycles.dat, and that the first column contains gnp * * in this section, the user can, of course, do any transformations * of the data * * !!! The program assumes, however, that all observations in * the matrix rawdat(ntot,nvar) are observations that can be used * to construct the filtered data!!! * * ******************************************************************** * * * THE USER DOES NOT HAVE TO CHANGE ANYTHING BELOW * * ******************************************************************** declare vect[series] rawser(nvar) declare vect[series] hpser(nvar) declare vect[series] trser(nvar) declare vect[series] bkser(nvar) declare vect[series] difser(nvar) declare vect[label] xstr(25) declare rectangle raw1(nvar,8) hp1(nvar,8) bk1(nvar,8) dif1(nvar,8) $ raw2(nvar,8) hp2(nvar,8) bk2(nvar,8) dif2(nvar,8) $ raw3(nvar,8) hp3(nvar,8) bk3(nvar,8) dif3(nvar,8) input xstr y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 $ x19 x20 x21 x22 x23 x24 x25 smpl nfirst nlast do i = 1,nvar set rawser(i) = 0. do j = nfirst,nlast compute rawser(i)(j) = rawdat(j+1-nfirst,i) end do j end do i *********************************************************************** * * DO THE PROCEDURES * *********************************************************************** PROCEDURE VARHAC start end *********************************************************************** *** *** DOS VERSION *** *** To use interactively (i.e., menues): *** *** @VARHAC start end *** # the dependent variables *** *** NOTE THAT YOU ANSWERS QUESTIONS BY TYPING THE ITEM NUMBER!!! *** *** *** To use noninteractively: *** *** @VARHAC(nointeractive,det= ,crit= ,maxlag= ) start end *** # the dependent variables *** *** where: *** *** det = "none" or "constant" -- deterministic component in the VAR *** crit = "AIC", "SIC" or "Kbar" -- model selection criteria *** (Kbar means use Kbar as the lag in each eq. in the VAR) *** maxlag = maxlag to consider for the model selection criteria *** *** *** In BOTH instances NOPRINT will suppress printed output *** *** *** The Spectral density estimator is saved to the symmetric matrix SDM *** *** References: *** Den Haan, Wouter J., and Andrew Levin, 1994: "Inferences from *** Parametric and Non-parametric Covariance Matrix *** Estimation Procedures" UCSD Discussion Paper, 94-30. *** *** *** EXAMPLES *** *** *** *** source(noecho) VARHAC.SRC *** *** @varhac(nointeractive,det=constant,crit=aic,maxlag=4) 59:1 93:4 *** #u1 u2 u3 u4 u5 *** *** *** or *** *** @varhac(nointeractive,noprint,det=none,crit=aic,maxlag=4) 59:1 93:4 *** #u1 u2 u3 u4 u5 *** *** *** *** !!! ALL VARIABLES IN THE VARHAC PROCEDURE ARE LOCAL. HOWEVER, *** THE VARHAC PROCEDURE ESTIMATES A VAR SYSTEM, AND THE PARAMETERS *** ESTIMATED BY THIS SYSTEM ARE NOT. *** *** ********************************************************************* DECLARE SYMM VCOVMAT SDM TYPE INTEGER start end OPTION CHOICE DET 1 NONE CONSTANT OPTION INTEGER MAXLAG 4 OPTION CHOICE CRIT 1 AIC SIC Kbar OPTION SWITCH INTERACTIVE 1 OPTION SWITCH PRINT 1 LOCAL INDEX varlist levsx ccmat laglevsx LOCAL INTEGER lag vars counter startl endl nprint num deter critmax LOCAL INTEGER Kbar Kmax ii jj kk i j k nvars msc LOCAL RECT lev laglev f11 F FP LOCAL VECT[series] DEPVAR RESIDVEC RESIDS LOCAL RECT cmat AIC SIC U Y YLAGS LOCAL SYMM VCOV VCORR SDMAT LOCAL VECT bb bb2 STDEV RSSMAT LOCAL SERIES AICVEC SICVEC LOCAL VECT AICMINENT SICMINENT AICMIN SICMIN LOCAL VECT[rect] outmatrix LOCAL VECT[labels] vlabel LOCAL VECT[rect] VMAMAT COEFFMAT LOCAL VECT[INDEX] laglevsx2 LOCAL VECT[VECT] COEFFVEC ENTER(varying) varlist COMPUTE vars = %rows(varlist) INQUIRE(regressorlist) startl>>start endl>>end #varlist DIM vlabel(vars) DO ii = 1,vars COMPUTE vlabel(ii) = %label([series] varlist(ii)) END DO IF PRINT.eq.1 { DISP; DISP DISP '*******************************************************' DISP '* The VARHAC Estimator of the Spectral Density Matrix *' DISP '* by Norman Morin *' DISP '*******************************************************' DISP } DIMENSION depvar(vars) DIMENSION residvec(vars) DIMENSION coeffvec(vars) DIMENSION levsx(0) DIMENSION laglevsx(0) DIMENSION ccmat(0) IF INTERACTIVE.eq.0 { COMPUTE Kbar = maxlag COMPUTE deter = det COMPUT MSC = crit } ELSE { MENU 'Deterministic components' CHOICE 'None' COMPUTE deter = 1 CHOICE 'Constant' COMPUTE deter = 2 END MENU DISPLAY(store=qstring) $ 'Input the maximum lag to consider for AIC/Schwarz' QUERY(prompt=qstring,status=check) Kbar MENU 'Model Selection Criterion' CHOICE 'AIC' COMPUTE MSC = 1 CHOICE 'Schwarz' COMPUTE MSC = 2 CHOICE 'The maximum lag previously input' COMPUTE MSC = 3 END MENU } DO i = 1,vars ENTER(varying) levsx # levsx varlist(i) ENTER(varying) laglevsx # laglevsx varlist(i){1 to Kbar} END DO MAKE(trans) lev startl endl # levsx MAKE(trans) laglev startl endl # laglevsx DO ii = 1,vars SET depvar(ii) startl endl = lev(ii,t-startl+1) LABLES DEPVAR(ii) #VLABEL(ii) END DO IF MSC.lt.3 { ********************************************* IF (deter.eq.1) { DIM aic(vars,Kbar) COMPUTE aic = %CONST(0.) DIM sic(vars,Kbar) COMPUTE sic = %CONST(0.) DO lagnum = 1,Kbar DIM RESIDVEC(vars) SYSTEM 1 to vars LAGS 1 to lagnum VARIABLES depvar END(SYSTEM) ESTIMATE(noftest,noprint) startl+Kbar endl RESIDVEC(1) DO jj = 1,vars set rr1 startl+Kbar endl = RESIDVEC(jj)**2 acc rr1 startl+Kbar endl rr2 COMPUTE aic(jj,lagnum) = $ log(rr2(endl)/%nobs) + 2.*%nreg/%nobs COMPUTE sic(jj,lagnum) = $ log(rr2(endl)/%nobs) + log(%nobs)*(1.*%nreg/%nobs) END DO END DO } ELSE { DIM aic(vars,Kbar) COMPUTE aic = %CONST(0.) DIM sic(vars,Kbar) COMPUTE sic = %CONST(0.) DO lagnum = 1,Kbar DIM RESIDVEC(vars) SYSTEM 1 to vars LAGS 1 to lagnum VARIABLES depvar DETERMINISTIC CONSTANT END(SYSTEM) ESTIMATE(noftest,noprint) startl+Kbar endl RESIDVEC(1) DO jj = 1,vars set rr1 startl+Kbar endl = RESIDVEC(jj)**2 acc rr1 startl+Kbar endl rr2 COMPUTE aic(jj,lagnum) = $ log(rr2(endl)/%nobs) + 2.*%nreg/%nobs COMPUTE sic(jj,lagnum) = $ log(rr2(endl)/%nobs) + log(%nobs)*(1.*%nreg/%nobs) END DO END DO } SET AICVEC 1 Kbar = 0.; SET SICVEC 1 Kbar = 0. DIM AICMIN(vars) SICMIN(vars) AICMINENT(vars) SICMINENT(vars) DO jj = 1,vars DO ii = 1,Kbar COMPUTE aicvec(ii) = aic(jj,ii) COMPUTE sicvec(ii) = sic(jj,ii) END DO EXTREMUM(noprint) aicvec 1 Kbar COMPUTE aicminent(jj) = %minent COMPUTE aicmin(jj) = %minimum EXTREMUM(noprint) sicvec 1 Kbar COMPUTE sicminent(jj) = %minent COMPUTE sicmin(jj) = %minimum END DO IF MSC.eq.1 { COMPUTE critmax = FIX(AICMINENT(1)) DO jj = 1,vars IF critmax.lt.(FIX(AICMINENT(jj))) COMPUTE critmax = FIX(AICMINENT(jj)) END DO } ELSE { COMPUTE critmax = FIX(SICMINENT(1)) DO jj = 1,vars IF critmax.lt.(FIX(SICMINENT(jj))) COMPUTE critmax = FIX(SICMINENT(jj)) END DO } IF PRINT.eq.1 { DISP ' The maximum lag length considered is' Kbar DISP WRITE(noskip,format='(7X8A9)') VLABEL DISP ' ' DISP 'SIC' WRITE(noskip,format='(8F9.0)') SICMINENT DISP 'AIC' WRITE(noskip,format='(8F9.0)') AICMINENT DISP IF MSC.eq.1 DISP ' The AIC is being used as the model selection criterion' ELSE DISP ' The SIC is being used as the model selection criterion' DISP } } **************************** ELSE { IF PRINT.eq.1 { DISP ' Using' Kbar 'lags in each equation in the VAR' DISP } DIM AICMIN(vars) SICMIN(vars) AICMINENT(vars) SICMINENT(vars) COMPUTE AICMINENT = %const(kbar) COMPUTE SICMINENT = %const(kbar) COMPUTE critmax = kbar } ***************************************************************** DIM laglevsx2(vars) RESIDS(vars) COEFFVEC(vars) DO jj = 1,vars DIM COEFFVEC(jj)(FIX(AICMINENT(jj))) DIM LAGLEVSX2(jj)(0) DO i = 1,vars IF (MSC.eq.1) ENTER(varying) laglevsx2(jj) # laglevsx2(jj) varlist(i){1 to FIX(AICMINENT(jj))} ELSE ENTER(varying) laglevsx2(jj) # laglevsx2(jj) varlist(i){1 to FIX(SICMINENT(jj))} END DO IF (DETER.eq.1) { LINREG(noprint) depvar(jj) startl+Kbar endl RESIDS(jj) #laglevsx2(jj) COMPUTE COEFFVEC(jj) = %beta } ELSE { LINREG(noprint) depvar(jj) startl+Kbar endl RESIDS(jj) #laglevsx2(jj) constant COMPUTE COEFFVEC(jj) = %beta } END DO DIM U(vars,endl-startl-Kbar+1) Y(vars,endl-startl-Kbar+1) DO jj = 1,vars DO ii = startl+Kbar,endl COMPUTE U(jj,ii-(startl+Kbar)+1) = RESIDS(jj)(ii) COMPUTE Y(jj,ii-(startl+Kbar)+1) = DEPVAR(jj)(ii) END DO END DO COMPUTE nin = 1./(endl-(startl+Kbar)+1) COMPUTE VCOVMAT = nin*(U*tr(U)) DIM COEFFMAT(critmax) DO jj = 1,critmax DIM COEFFMAT(jj)(vars,vars) COMPUTE COEFFMAT(jj) = %CONST(0.) END DO IF MSC.eq.1 { DO ii = 1,vars ; ******** EQ # DO kk = 1,FIX(AICMINENT(ii)) ; ******** # LAGS DO jj = 1,vars ; ******** # of VARIABLE in iith EQ COMPUTE COEFFMAT(kk)(ii,jj) = $ coeffvec(ii)(kk+(jj-1)*FIX(AICMINENT(ii)) ) END DO END DO END DO } ELSE { DO ii = 1,vars ; ******** EQ # DO kk = 1,FIX(SICMINENT(ii)) ; ******** # LAGS DO jj = 1,vars ; ******** # of VARIABLE in iith EQ COMPUTE COEFFMAT(kk)(ii,jj) = $ coeffvec(ii)(kk+(jj-1)*FIX(SICMINENT(ii)) ) END DO END DO END DO } COMPUTE SUMMAT = %IDENTITY(vars) COMPUTE SUMMAT2 = %IDENTITY(vars) DO jj = 1,critmax COMPUTE SUMMAT = SUMMAT - COEFFMAT(jj) END DO COMPUTE SDM = INV(SUMMAT)*VCOVMAT*TR(INV(SUMMAT)) IF PRINT.eq.1 { DISP ' ' DISP 'The Estimated Spectral Density Matrix at the Zero Frequency' DISP ' ' WRITE(noskip,format='(7X8A14)') VLABEL WRITE(noskip,format='(8g14.5)') SDM } END * * program to calculate the standard deviation * and its standard error for a series * * * you have to use the VARHAC procedure as well * * PROCEDURE TSTATSD Y1 BDATA EDATA SD TSTAT1 TSTAT2 type SERIES Y1 type INTEGER BDATA EDATA TYPE REAL *SD TYPE REAL *TSTAT1 TYPE REAL *TSTAT2 LOCAL SERIES COVARS COVS smpl bdata edata statistics(noprint) y1 / compute mean1 = %mean compute sd = sqrt(%variance) set res = (y1(t)-mean1)**2 - sd**2. * * * DEN HAAN and LEVIN "VARHAC" PROCEDURE * * @VARHAC(nointeractive,noprint,det=constant,crit=AIC,maxlag=4) bdata edata #res compute nt = edata-bdata+1 compute nn = nt**0.5 compute sdm(1,1) = sdm(1,1)/(4.*sd*sd) compute tstat1 = nn*sd/(sdm(1,1)**0.5) * * * NEWEY & WEST "OPTIMAL BANDWIDTH" PROCEDURE * * * * * first step calculate the optimal bandwidth * * compute dnn = 4.*((nt/100.)**(2./9.)) compute [integer] inn = fix(dnn) correlate(number=inn,covariances,noprint) res bdata edata covars compute spec0 = covars(1) compute spec1 = 0. do ii = 1,inn compute spec0 = spec0 + 2.*covars(1+ii) compute spec1 = spec1 + 2.*covars(1+ii)*ii end do ii compute zalpha = (spec1/spec0)**2. compute ststar = 1.1447*(zalpha*nt)**(1./3.) compute [integer] inn = fix(ststar) compute ststar = 1.*inn + 1. * * * second step calculate the zero-frequency optimal bandwidth * * correlate(number=inn,covariances,noprint) res bdata edata covs compute spec0 = covs(1) do ii = 1,inn compute weight = (1.-ii/ststar) compute spec0 = spec0 + 2.*weight*covs(1+ii) end do ii compute spec0 = spec0/(4.*sd*sd) compute tstat2 = nn*sd/(spec0**0.5) END * * program to calculate the standard error for * the correlation coefficient between Y1(t) and y2(t) * * * you have to use the VARHAC procedure as well * * PROCEDURE TSTATCOR Y1 Y2 BDATA EDATA CORREL TSTAT1 TSTAT2 TYPE SERIES Y1 Y2 TYPE INTEGER BDATA EDATA TYPE REAL *CORREL TYPE REAL *TSTAT1 TYPE REAL *TSTAT2 LOCAL SERIES COVARS COVS LOCAL RECT D DIMENSION D(2,2) smpl bdata edata statistics(noprint) y1 / compute mean1 = %mean compute sd1 = sqrt(%variance) statistics(noprint) y2 / compute mean2 = %mean compute sd2 = sqrt(%variance) compute estimate = %corr(y1,y2) set res1 = (y1(t)-mean1)**2.*(sd2/sd1)**2. - (y2(t)-mean2)**2. set res2 = (y1(t)-mean1)**2.*(sd2/sd1)*estimate - $ (y2(t)-mean2)*(y1(t)-mean1) compute d(1,1) = 2.*sd1*sd2 compute d(1,2) = 0.0 compute d(2,1) = estimate*sd1*sd1 compute d(2,2) = sd1*sd2 * * * DEN HAAN and LEVIN "VARHAC" PROCEDURE * * @VARHAC(nointer,noprint,det=constant,crit=AIC,maxlag=4) bdata edata #res1 res2 compute nt = edata-bdata+1 compute nn = nt**0.5 compute vvv = inv(d)*sdm*inv(tr(d)) compute tstat1 = nn*estimate/(vvv(2,2)**0.5) * * * NEWEY & WEST "OPTIMAL BANDWIDTH" PROCEDURE * * * * * first step calculate the optimal bandwidth, note that * we only use the persistence in the second residual * to determine the optimal bandwidth. * compute dnn = 4.*((nt/100.)**(2./9.)) compute [integer] inn = fix(dnn) correlate(number=inn,covariances,noprint) res2 bdata edata covars compute spec0 = covars(1) compute spec1 = 0. do ii = 1,inn compute spec0 = spec0 + 2.*covars(1+ii) compute spec1 = spec1 + 2.*covars(1+ii)*ii end do ii compute zalpha = (spec1/spec0)**2. compute ststar = 1.1447*(zalpha*nt)**(1./3.) compute [integer] inn = fix(ststar) compute ststar = 1.*inn + 1 * * * second step calculate the zero-frequency optimal bandwidth * * compute sdm(1,1) = 0. compute sdm(1,2) = 0. compute sdm(2,1) = 0. compute sdm(2,2) = 0. do i = bdata, edata compute sdm(1,1) = sdm(1,1) + res1(i)*res1(i) compute sdm(1,2) = sdm(1,2) + res1(i)*res2(i) compute sdm(2,1) = sdm(2,1) + res1(i)*res2(i) compute sdm(2,2) = sdm(2,2) + res2(i)*res2(i) end do i do ii = 1,inn compute weight = (1.-ii/ststar) do i = bdata+ii, edata compute sdm(1,1) = sdm(1,1) + res1(i)*res1(i-ii)*2.*weight compute sdm(1,2) = sdm(1,2) + res1(i)*res2(i-ii)*2.*weight compute sdm(2,1) = sdm(2,1) + res2(i)*res1(i-ii)*2.*weight compute sdm(2,2) = sdm(2,2) + res2(i)*res2(i-ii)*2.*weight end do i end do ii compute sdm(1,1) = sdm(1,1)/nt compute sdm(1,2) = sdm(1,2)/nt compute sdm(2,2) = sdm(2,2)/nt compute sdm(2,1) = sdm(2,1)/nt compute vvv = inv(d)*sdm*inv(tr(d)) compute tstat2 = nn*estimate/(vvv(2,2)**0.5) compute correl = estimate END * * program to calculate the standard error for * the ratio of sd(y1(t)) divided by sd(y2(t)) * * * * you have to use the VARHAC procedure as well * * PROCEDURE TSTATRAT Y1 Y2 BDATA EDATA ESTIMATE TSTAT1 TSTAT2 TYPE SERIES Y1 Y2 TYPE INTEGER BDATA EDATA TYPE REAL *ESTIMATE TYPE REAL *TSTAT1 TYPE REAL *TSTAT2 LOCAL SERIES COVARS COVS smpl bdata edata statistics(noprint) y1 / compute mean1 = %mean compute sd1 = sqrt(%variance) statistics(noprint) y2 / compute mean2 = %mean compute sd2 = sqrt(%variance) compute estimate = sd1/sd2 set res = (y1(t)-mean1)**2.*(sd2/sd1)**2. - (y2(t)-mean2)**2. * * * DEN HAAN and LEVIN "VARHAC" PROCEDURE * * @VARHAC(nointeractive,noprint,det=constant,crit=AIC,maxlag=4) bdata edata #res compute nt = edata-bdata+1 compute nn = nt**0.5 compute se = sdm(1,1)/(4.*(sd1*sd2)**2.) compute tstat1 = nn*estimate/(se**0.5) * * * * NEWEY & WEST "OPTIMAL BANDWIDTH" PROCEDURE * * * * * first step calculate the optimal bandwidth * * compute dnn = 4.*((nt/100.)**(2./9.)) compute [integer] inn = fix(dnn) correlate(number=inn,covariances,noprint) res bdata edata covars compute spec0 = covars(1) compute spec1 = 0. do ii = 1,inn compute spec0 = spec0 + 2.*covars(1+ii) compute spec1 = spec1 + 2.*covars(1+ii)*ii end do ii compute zalpha = (spec1/spec0)**2. compute ststar = 1.1447*(zalpha*nt)**(1./3.) compute [integer] inn = fix(ststar) compute ststar = 1.*inn + 1. * * * second step calculate the zero-frequency optimal bandwidth * * correlate(number=inn,covariances,noprint) res bdata edata covs compute spec0 = covs(1) do ii = 1,inn compute weight = (1.-ii/ststar) compute spec0 = spec0 + 2.*weight*covs(1+ii) end do ii compute se = spec0/(4.*(sd1*sd2)**2.) compute tstat2 = nn*estimate/(se**0.5) END * * HPFILTER.SRC * Executes a Hodrick-Prescott Filter (Hodrick, R. & Prescott, E., * "Post-War U.S. Business Cycles: An Empirical Investigation", * Carnegie-Mellon working paper, 1980. * * Syntax: * * @HPFILTER series start end growth component * * Options: * LAMBDA=weight on squared 2nd difference of growth component [1600] * * HPFILTER computes the estimated growth component using the DP algorithm * to solve the problem: * * Minimize sum over t of * * (*) (y(t)-g(t))**2 + lambda * (g(t)-2g(t-1)+g(t-2))**2 * * where g is the growth component of y. * * X(t)=(g(t),g(t-1),1)' * u(t)=g(t)-2g(t-1)+g(t-2) * * X(t) = A X(t-1) + B u(t), where * * 2 -1 0 1 * A = 1 0 0 B = 0 * 0 0 1 0 * * Using these definitions, the minimand (*) can be written * * (c(t)X(t))**2 + lambda * u(t)**2 * c(t)=(1,0,-y(t)) * * With cost-to-go functional X(t)'R(t)X(t), the solution is of the form * u(t)=TEMP1(t)*X(t-1), and X(t)=(A+B*TEMP1(t))*X(t-1) * PROCEDURE HPFILTER SERIES START END OUTSER TYPE SERIES SERIES TYPE INTEGER START END TYPE SERIES *OUTSER * OPTION REAL LAMBDA 1600 * LOCAL INTEGER STARTL ENDL * LOCAL RECT A B C LOCAL SYMM R RR LOCAL RECT D LOCAL RECT DMATS LOCAL RECT TEMP1 LOCAL VECT RB X LOCAL INTEGER TOTAL I J * IF .NOT.(%DEFINED(SERIES).AND.%DEFINED(OUTSER)) { DISPLAY 'Syntax: @HPFILTER series start end growth component' RETURN } INQUIRE(SERIES=SERIES) STARTL>>START ENDL>>END COMPUTE TOTAL=ENDL-STARTL DIM A(3,3) B(3,1) C(1,3) R(3,3) X(3) DIM DMATS(9,TOTAL) COMPUTE A=||2,-1,0|1,0,0|0,0,1|| COMPUTE B=||1|0|0|| COMPUTE C=||1,0,0|| COMPUTE R=%CONST(0.0) * * DP transition matrices are set up. Loop over periods until 2 above * beginning of data. * DO I=ENDL,STARTL+2,-1 COMPUTE C(1,3)=-SERIES(I) OVERLAY DMATS(1,I-STARTL) WITH D(3,3) COMPUTE R=R+TR(C)*C COMPUTE SCALAR=LAMBDA+%QFORM(R,B) COMPUTE TEMP1=(-1.0/SCALAR)*TR(B)*R*A COMPUTE D=A+B*TEMP1 COMPUTE R=%MQFORM(R,D)+LAMBDA*TR(TEMP1)*TEMP1 END DO I * * Compute minimizer for the initial two periods. Since s(1) and s(0) * are free initial parameters, u(1) and u(2) can be set equal to zero. * Thus, we only need to deal with the (y(t)-s(t))**2 terms. * COMPUTE C(1,3)=-SERIES(STARTL+1) COMPUTE R=R+TR(C)*C COMPUTE C=||0.0,1.0,-SERIES(STARTL)|| COMPUTE R=R+TR(C)*C * OVERLAY R(1,1) WITH RR(2,2) OVERLAY R(3,1) WITH RB(2) * SET OUTSER STARTL ENDL = 0.0 COMPUTE RB=INV(RR)*RB COMPUTE X(1)=OUTSER(STARTL+1)=-RB(1) COMPUTE X(2)=OUTSER(STARTL)=-RB(2) COMPUTE X(3)=1.0 * DO I=STARTL+2,ENDL OVERLAY DMATS(1,I-STARTL) WITH D(3,3) COMPUTE X=D*X COMPUTE OUTSER(I)=X(1) END DO I END procedure filter1 xin ss ee number kkk xout type real number type integer ss type integer ee type integer kkk type series xin type series *xout local integer start local integer end local vector b local vector a local vector temp dimension b(1000) a(1000) temp(300000) do i = 1,1000 compute b(i) = 0. compute a(i) = 0. end do i compute wbar = 2.*%pi/number compute [integer] start = ss+kkk compute [integer] end = ee-kkk compute b0 = wbar/%pi do i = 1,kkk compute b(i) = sin(wbar*i)/(%pi*i) end do i * * note that the correction is as in B&K(1994) not B&K(1993) * compute sumb = %sum(b) compute sumb = b0+2.*sumb compute theta = (1.-sumb)/(2.*kkk+1) compute a0 = b0 + theta do i = 1,kkk compute a(i) = b(i) + theta end do i do time = start, end compute xxx = xin(time)*a0 do i = 1,kkk compute xxx = xxx + a(i)*( xin(time+i)+xin(time-i) ) end do i compute temp(time) = xin(time) - xxx end do time smpl start end set xout = temp(t) end * * * FILTER THE SERIES * * do iser = 1,nvar if ihp == 1 { smpl nfirst nlast @hpfilter rawser(iser) nfirst nlast trser(iser) set hpser(iser) = rawser(iser) - trser(iser) } if ibk == 1 { smpl nfirst nlast @filter1 rawser(iser) nfirst nlast bkper bklost bkser(iser) } { if idif == 1 smpl nfirst+1 nlast set difser(iser) = rawser(iser)(t) - rawser(iser)(t-1) } end do iser * * * CALCULATE STATISTICS IF IRAW = 1 * * if iraw == 1 { smpl start end @tstatsd rawser(1) start end raw1(1,1) raw2(1,1) raw3(1,1) do ilag = 1,3 compute [integer] ss = start+ilag smpl ss end set t1 = rawser(1)(t-ilag) @tstatcor rawser(1) t1 ss end raw1(1,5-ilag) raw2(1,5-ilag) raw3(1,5-ilag) compute [integer] ee = end - ilag smpl start ee set t1 = rawser(1)(t+ilag) @tstatcor rawser(1) t1 start ee raw1(1,5+ilag) raw2(1,5+ilag) raw3(1,5+ilag) end do ilag compute raw1(1,5) = 1. do iv = 2,nvar smpl start end @tstatrat rawser(iv) rawser(1) start end raw1(iv,1) raw2(iv,1) raw3(iv,1) @tstatcor rawser(iv) rawser(1) start end raw1(iv,5) raw2(iv,5) raw3(iv,5) do ilag = 1,3 compute [integer] ss = start+ilag smpl ss end set t1 = rawser(iv)(t-ilag) @tstatcor rawser(1) t1 ss end raw1(iv,5-ilag) raw2(iv,5-ilag) raw3(iv,5-ilag) compute [integer] ee = end - ilag smpl start ee set t1 = rawser(iv)(t+ilag) @tstatcor rawser(1) t1 start ee raw1(iv,5+ilag) raw2(iv,5+ilag) raw3(iv,5+ilag) end do ilag end do iv * * * transform t-statistics back into standard errors * * do i = 1,nvar do j = 1,8 compute raw2(i,j) = raw1(i,j)/raw2(i,j) compute raw3(i,j) = raw1(i,j)/raw3(i,j) end do j end do i DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' BUSINESS CYCLE STATISTICS FOR RAW DATA ' DISPLAY(UNIT=COPY) ' FROM ' %datelabel(START) ' TO ' %datelabel(END) DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### raw1(i,1) $ @+1 ##.#### raw1(i,2) @+1 ##.#### raw1(i,3) $ @+1 ##.#### raw1(i,4) @+1 ##.#### raw1(i,5) $ @+1 ##.#### raw1(i,6) @+1 ##.#### raw1(i,7) $ @+1 ##.#### raw1(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' VARHAC STANDARD ERRORS ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### raw2(i,1) $ @+1 ##.#### raw2(i,2) @+1 ##.#### raw2(i,3) $ @+1 ##.#### raw2(i,4) @+1 ##.#### raw2(i,5) $ @+1 ##.#### raw2(i,6) @+1 ##.#### raw2(i,7) $ @+1 ##.#### raw2(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' NEWEY & WEST OPTIMAL-BANDWIDTH STANDARD ERRORS ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### raw3(i,1) $ @+1 ##.#### raw3(i,2) @+1 ##.#### raw3(i,3) $ @+1 ##.#### raw3(i,4) @+1 ##.#### raw3(i,5) $ @+1 ##.#### raw3(i,6) @+1 ##.#### raw3(i,7) $ @+1 ##.#### raw3(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' } * * * CALCULATE STATISTICS IF IHP = 1 * * if ihp == 1 { smpl start end @tstatsd hpser(1) start end hp1(1,1) hp2(1,1) hp3(1,1) do ilag = 1,3 compute [integer] ss = start+ilag smpl ss end set t1 = hpser(1)(t-ilag) @tstatcor hpser(1) t1 ss end hp1(1,5-ilag) hp2(1,5-ilag) hp3(1,5-ilag) compute [integer] ee = end - ilag smpl start ee set t1 = hpser(1)(t+ilag) @tstatcor hpser(1) t1 start ee hp1(1,5+ilag) hp2(1,5+ilag) hp3(1,5+ilag) end do ilag compute hp1(1,5) = 1. do iv = 2,nvar smpl start end @tstatrat hpser(iv) hpser(1) start end hp1(iv,1) hp2(iv,1) hp3(iv,1) @tstatcor hpser(iv) hpser(1) start end hp1(iv,5) hp2(iv,5) hp3(iv,5) do ilag = 1,3 compute [integer] ss = start+ilag smpl ss end set t1 = hpser(iv)(t-ilag) @tstatcor hpser(1) t1 ss end hp1(iv,5-ilag) hp2(iv,5-ilag) hp3(iv,5-ilag) compute [integer] ee = end - ilag smpl start ee set t1 = hpser(iv)(t+ilag) @tstatcor hpser(1) t1 start ee hp1(iv,5+ilag) hp2(iv,5+ilag) hp3(iv,5+ilag) end do ilag end do iv * * * transform t-statistics back into standard errors * * do i = 1,nvar do j = 1,8 compute hp2(i,j) = hp1(i,j)/hp2(i,j) compute hp3(i,j) = hp1(i,j)/hp3(i,j) end do j end do i DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' BUSINESS CYCLE STATISTICS FOR HP FILTERED DATA ' DISPLAY(UNIT=COPY) ' FROM ' %datelabel(START) ' TO ' %datelabel(END) DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### hp1(i,1) $ @+1 ##.#### hp1(i,2) @+1 ##.#### hp1(i,3) $ @+1 ##.#### hp1(i,4) @+1 ##.#### hp1(i,5) $ @+1 ##.#### hp1(i,6) @+1 ##.#### hp1(i,7) $ @+1 ##.#### hp1(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' VARHAC STANDARD ERRORS ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### hp2(i,1) $ @+1 ##.#### hp2(i,2) @+1 ##.#### hp2(i,3) $ @+1 ##.#### hp2(i,4) @+1 ##.#### hp2(i,5) $ @+1 ##.#### hp2(i,6) @+1 ##.#### hp2(i,7) $ @+1 ##.#### hp2(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' NEWEY & WEST OPTIMAL-BANDWIDTH STANDARD ERRORS ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### hp3(i,1) $ @+1 ##.#### hp3(i,2) @+1 ##.#### hp3(i,3) $ @+1 ##.#### hp3(i,4) @+1 ##.#### hp3(i,5) $ @+1 ##.#### hp3(i,6) @+1 ##.#### hp3(i,7) $ @+1 ##.#### hp3(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' } * * * CALCULATE STATISTICS IF IBK = 1 * * if ibk == 1 { compute [integer] xstart = start if start.lt.nfirst+bklost { compute [integer] xstart = nfirst + bklost } compute [integer] xend = end if end.gt.nlast-bklost { compute [integer] xend = nlast - bklost } @tstatsd bkser(1) xstart xend bk1(1,1) bk2(1,1) bk3(1,1) do ilag = 1,3 compute [integer] ss = xstart+ilag smpl ss xend set t1 = bkser(1)(t-ilag) @tstatcor bkser(1) t1 ss xend bk1(1,5-ilag) bk2(1,5-ilag) bk3(1,5-ilag) compute [integer] ee = xend - ilag smpl xstart ee set t1 = bkser(1)(t+ilag) @tstatcor bkser(1) t1 xstart ee bk1(1,5+ilag) bk2(1,5+ilag) bk3(1,5+ilag) end do ilag compute bk1(1,5) = 1. do iv = 2,nvar smpl xstart xend @tstatrat bkser(iv) bkser(1) xstart xend bk1(iv,1) bk2(iv,1) bk3(iv,1) @tstatcor bkser(iv) bkser(1) xstart xend bk1(iv,5) bk2(iv,5) bk3(iv,5) do ilag = 1,3 compute [integer] ss = xstart+ilag smpl ss xend set t1 = bkser(iv)(t-ilag) @tstatcor bkser(1) t1 ss xend bk1(iv,5-ilag) bk2(iv,5-ilag) bk3(iv,5-ilag) compute [integer] ee = xend - ilag smpl xstart ee set t1 = bkser(iv)(t+ilag) @tstatcor bkser(1) t1 xstart ee bk1(iv,5+ilag) bk2(iv,5+ilag) bk3(iv,5+ilag) end do ilag end do iv * * * transform t-statistics back into standard errors * * do i = 1,nvar do j = 1,8 compute bk2(i,j) = bk1(i,j)/bk2(i,j) compute bk3(i,j) = bk1(i,j)/bk3(i,j) end do j end do i DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' BUSINESS CYCLE STATISTICS FOR BK FILTERED DATA ' DISPLAY(UNIT=COPY) ' FROM ' %datelabel(XSTART) ' TO ' %datelabel(XEND) DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### bk1(i,1) $ @+1 ##.#### bk1(i,2) @+1 ##.#### bk1(i,3) $ @+1 ##.#### bk1(i,4) @+1 ##.#### bk1(i,5) $ @+1 ##.#### bk1(i,6) @+1 ##.#### bk1(i,7) $ @+1 ##.#### bk1(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' VARHAC STANDARD ERRORS FOR BK FILTERED DATA ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### bk2(i,1) $ @+1 ##.#### bk2(i,2) @+1 ##.#### bk2(i,3) $ @+1 ##.#### bk2(i,4) @+1 ##.#### bk2(i,5) $ @+1 ##.#### bk2(i,6) @+1 ##.#### bk2(i,7) $ @+1 ##.#### bk2(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' NEWEY & WEST OPTIMAL-BANDWIDTH STANDARD ERRORS ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### bk3(i,1) $ @+1 ##.#### bk3(i,2) @+1 ##.#### bk3(i,3) $ @+1 ##.#### bk3(i,4) @+1 ##.#### bk3(i,5) $ @+1 ##.#### bk3(i,6) @+1 ##.#### bk3(i,7) $ @+1 ##.#### bk3(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' } * * * CALCULATE STATISTICS IF IDIF = 1 * * if idif == 1 { compute [integer] xstart = start if start.lt.nfirst+1 { compute [integer] xstart = nfirst + 1 } compute [integer] xend = end @tstatsd difser(1) xstart xend dif1(1,1) dif2(1,1) dif3(1,1) do ilag = 1,3 compute [integer] ss = xstart+ilag smpl ss xend set t1 = difser(1)(t-ilag) @tstatcor difser(1) t1 ss xend dif1(1,5-ilag) dif2(1,5-ilag) dif3(1,5-ilag) compute [integer] ee = xend - ilag smpl xstart ee set t1 = difser(1)(t+ilag) @tstatcor difser(1) t1 xstart ee dif1(1,5+ilag) dif2(1,5+ilag) dif3(1,5+ilag) end do ilag compute dif1(1,5) = 1. do iv = 2,nvar smpl xstart xend @tstatrat difser(iv) difser(1) xstart xend dif1(iv,1) dif2(iv,1) dif3(iv,1) @tstatcor difser(iv) difser(1) xstart xend dif1(iv,5) dif2(iv,5) dif3(iv,5) do ilag = 1,3 compute [integer] ss = xstart+ilag smpl ss xend set t1 = difser(iv)(t-ilag) @tstatcor difser(1) t1 ss xend dif1(iv,5-ilag) dif2(iv,5-ilag) dif3(iv,5-ilag) compute [integer] ee = xend - ilag smpl xstart ee set t1 = difser(iv)(t+ilag) @tstatcor difser(1) t1 xstart ee dif1(iv,5+ilag) dif2(iv,5+ilag) dif3(iv,5+ilag) end do ilag end do iv * * * transform t-statistics back into standard errors * * do i = 1,nvar do j = 1,8 compute dif2(i,j) = dif1(i,j)/dif2(i,j) compute dif3(i,j) = dif1(i,j)/dif3(i,j) end do j end do i DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' BUSINESS CYCLE STATISTICS FIRST-DIFFERENCED DATA ' DISPLAY(UNIT=COPY) ' FROM ' %datelabel(XSTART) ' TO ' %datelabel(XEND) DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### DIF1(i,1) $ @+1 ##.#### DIF1(i,2) @+1 ##.#### DIF1(i,3) $ @+1 ##.#### DIF1(i,4) @+1 ##.#### DIF1(i,5) $ @+1 ##.#### DIF1(i,6) @+1 ##.#### DIF1(i,7) $ @+1 ##.#### DIF1(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' VARHAC STANDARD ERRORS FOR FIRST-DIFFERENCED DATA ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### DIF2(i,1) $ @+1 ##.#### DIF2(i,2) @+1 ##.#### DIF2(i,3) $ @+1 ##.#### DIF2(i,4) @+1 ##.#### DIF2(i,5) $ @+1 ##.#### DIF2(i,6) @+1 ##.#### DIF2(i,7) $ @+1 ##.#### DIF2(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' NEWEY & WEST OPTIMAL-BANDWIDTH STANDARD ERRORS ' DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' SD(X)/SD(Y)' $ ' CROSS-CORRELATION OF OUTPUT WITH X(T+J)' DISPLAY(UNIT=COPY) ' J=-3 J=-2 J=-1 '$ 'J=0 J=1 J=2 J=3' DO I=1,nvar DISPLAY(UNIT=COPY) XSTR(I) @5 ##.###### DIF3(i,1) $ @+1 ##.#### DIF3(i,2) @+1 ##.#### DIF3(i,3) $ @+1 ##.#### DIF3(i,4) @+1 ##.#### DIF3(i,5) $ @+1 ##.#### DIF3(i,6) @+1 ##.#### DIF3(i,7) $ @+1 ##.#### DIF3(i,8) END DO I DISPLAY(UNIT=COPY) ' ' DISPLAY(UNIT=COPY) ' ' } end