? ----- PROGRAM \BOOKLIM\CH7R1.LIM June 1997 ? Began as c:\booklim\ch7r2.lim ? This LIMDEP program does all the analysis in Chapter 7 of ? A.C. Cameron and Pravin K. Trivedi (1998), ? REGRESSION ANALYSIS OF COUNT DATA, ? Econometric Society Monograph No.30, Cambridge University Press. ? For more details on the data see files racd7.1st ? A title such as Table 7.2 (col.3) means Table 7.2 column 3 ? The program runs on DOS version of Limdep 7.0 (limdep.exe dated 3/28/97) ? It may not be completely backward compatible with Limdep 6, ? but most of the program was originally written for Limdep 6. ? The program is not written as efficiently written as possible, ? partly because it began with Limdep 6. ? It is written in, hopefully, a transparent manner. ? Also in places intermediate output is given to enable checks. ? The data is the same data as originally used in ? J. Kennan, "The Duration of Contract strikes in U.S. Manufacturing", ? Journal of Econometrics, 1985, Vol. 28, pp.5-28. ? Dependent variable: STRIKES is number of strikes ? Regressor variabel: OUTPUT is level of economic activity ? There are 108 observations from 1968(1) to 1976(12) ? +++++++++ Open output file and read in the data ? OPEN; OUTPUT=c:\racd\racd7p1.out $ READ; NVAR=4; NOBS=108; file=c:\racd\racd7.asc; NAMES = STRIKES,OUTPUT,TIME,CONSTANT $ ? +++++++++ Descriptive Statistics and histogram TITLE; * RACD Chapter 7, p.231: Actual frequency distribution for STRIKES.$ HISTOGRAM; RHS=STRIKES; INT=18 $ TITLE; * RACD Table 7.1, p.231: Descriptive Statistics for Y and X. $ DSTAT; RHS = STRIKES,OUTPUT,TIME,CONSTANT $ ? --------- Time Series Plot of Data from 1968(1) to 1976(12) TITLE; * RACD FIGURE 7.1, p.232: Plot of STRIKES and rescaled OUTPUT $ CREATE; ACTSCALE =(3.75131/0.054557)*OUTPUT + 5.24074+0.254088 $ ? A check that rescaled so ACTSCALE has same mean and variance as STRIKES DSTAT; RHS = STRIKES,ACTSCALE $ PLOT; RHS = STRIKES,ACTSCALE $ ? +++++++++ Poisson PMLE for static model (page 232) ? --------- Poisson PMLE and possibly incorrect ML standard errors TITLE; * RACD Table 7.2, col.1, p.232: Poisson PMLE $ NAMELIST; XSTATIC = ONE,OUTPUT $ POISSON; LHS=STRIKES; RHS=XSTATIC $ CALC; LIST; PLOGL = LOGL; PKREG = KREG; PN = N $ ? Get Pearson residual CREATE; PPRED = EXP(DOT(XSTATIC,B)); PPREDSQ = PPRED*PPRED; PRES = STRIKES - PPRED; PEARSRES = (STRIKES-PPRED)/SQR(PPRED) $ ? --------- Get NB1 standard errors (method of McCullagh-Nelder) TITLE; * RACD TABLE 7.2, cols.2-3, p.232: Poisson PMLE with NB1 s.e.'s $ ? First, GLM Estimate of NB1 Overdispersion parameter. CALC; NOLIST; DEGOFFR = PN - PKREG $ MATRIX; LIST; PEARSON=XDOT(PEARSRES); ANB1MN = {1/DEGOFFR}*PEARSON $ ? Second, ML standard errors and t-statistics MATRIX; NOLIST; UNITY = INIT(PKREG,1,1); DIAGV = DIAG(VARB); PSEML = SQRT(DIAGV) | UNITY; PTML = ISQR(DIAGV) | B $ ? Third, correct using estimated scale parameter CALC; TAU = ANB1MN $ MATRIX; NOLIST; TAUSQRT = SQRT(TAU) $ CALC; INVTAUSQ = 1/TAUSQRT $ MATRIX; LIST; PSEMN = TAUSQRT*PSEML; PTMN = INVTAUSQ*PTML $ ? +++++++++ Get other standard errors of static Poisson regression ? Both NLS and Poisson PML with exponential conditional mean function ? First het consistent uses pds = 0 ? I tried autocorrelation and het consistent with pds=1 and pds=5 ? as there is autocorrelation to 5 lags ? But LIMDEP reported singularity. SAMPLE; 1-108 $ TITLE; * RACD TABLE 7.2: Poisson PML het-consistent standard errors $ NLSQ; LHS = STRIKES; LABELS = BONE,BOUTPUT; START = 1.6,3.1; FCN = EXP(BONE + BOUTPUT*OUTPUT); INST = ONE,OUTPUT; pds=0 $ ? --------- Plot residuals and predicted values PLOT; RHS = STRIKES,PPRED $ PLOT; RHS = PRES $ ? The following is used later CREATE; ZQ0PRED = PPRED $ ? +++++++++ Autocorrelation functions of Poisson residuals (page 233) TITLE; * RACD TABLE 7.3, col.1, p.233: Autocorrelations of STRIKES $ IDENTIFY; RHS = STRIKES; pds = 15 $ TITLE; * RACD TABLE 7.3 col.2, p.233: Autocorrelations Poisson raw resid $ IDENTIFY; RHS = PRES; pds = 15 $ TITLE; * RACD TABLE 7.3 col.3, p.233: Autocorrelations Poisson Pearson res $ IDENTIFY; RHS = PEARSRES; pds = 15 $ ? --------- Autocorrelations using nonstandardized residuals ? This gives correct t-statistics for Poisson raw residuals ? from lag 1 to lag 6. ? Compare to other t-stats which are ACF times square root T TITLE; * RACD TABLE 7.3, p.233: Correct t-stats for Poisson raw residuals $ SAMPLE; 1-108 $ CREATE; R = PRES; RL1=R[-1]; RL2=R[-2]; RL3=R[-3]; RL4=R[-4]; RL5=R[-5]; RL6=R[-6]; RRL1 = R*RL1; RRL1SQ = RRL1*RRL1; RRL2 = R*RL2; RRL2SQ = RRL2*RRL2; RRL3 = R*RL3; RRL3SQ = RRL3*RRL3; RRL4 = R*RL4; RRL4SQ = RRL4*RRL4; RRL5 = R*RL5; RRL5SQ = RRL5*RRL5; RRL6 = R*RL6; RRL6SQ = RRL6*RRL6 $ ? Next is the comparison group. It uses residuals from 7-108 TITLE; * RACD TABLES 7.3, col.4: ACF Poisson raw res to lag 6$ ? Multiply these by square root of t to get t-statistic IDENTIFY; RHS = PEARSRES; pds = 6 $ TITLE; * RACD TABLE 7.3, col.5: ACF of Poisson Pearson residuals to lag 6$ SAMPLE; 7-108 $ CALC; LIST; T1 = SUM(RRL1)/SQR(SUM(RRL1SQ)); T2 = SUM(RRL2)/SQR(SUM(RRL2SQ)); T3 = SUM(RRL3)/SQR(SUM(RRL3SQ)); T4 = SUM(RRL4)/SQR(SUM(RRL4SQ)); T5 = SUM(RRL5)/SQR(SUM(RRL5SQ)); T6 = SUM(RRL6)/SQR(SUM(RRL6SQ)); BP6 = T1*T1 + T2*T2 + T3*T3 + T4*T4 + T5*T5 + T6*T6 $ ? +++++++++ Transform data for dynamic regressions SAMPLE; 1-108 $ ? In particular create yslag1 = max(c,ylag1) ? From now on use Y,Z terminology. CREATE; Y = STRIKES; Z = OUTPUT $ ? Next variable changes Y to 0.5 if Y = 0 CREATE; IF (Y = 0) YSTAR = 0.5; (else) YSTAR = Y $ CREATE; LNYSTAR = LOG(YSTAR) $ ? Now create lags, up to lag 3 for Y CREATE; YLAG1 = Y[-1]; LNYSLAG1 = LNYSTAR[-1]; YLAG2 = Y[-2]; LNYSLAG2 = LNYSTAR[-2]; YLAG3 = Y[-3]; LNYSLAG3 = LNYSTAR[-3] $ ? +++++++++ Dynamic regressions (page 247) ? --------- Initial analysis on how many lags are needed ? Use Zeger-Qaqish model with c = 0.5 ?TITLE; * RACD TABLE 7.4, col.1, p.247: Zeger-Qaqish no lag $ SAMPLE; 1-108 $ POISSON; LHS=STRIKES; RHS=XSTATIC; RES=RESID $ CREATE; PEARSRES = RESID/SQR(STRIKES-RESID) $ IDENTIFY; RHS = PEARSRES; pds = 6 $ TITLE; * RACD TABLE 7.4, col.2, p.247: Zeger-Qaqish 1 lag $ SAMPLE; 2-108 $ POISSON; LHS=STRIKES; RHS=XSTATIC,LNYSLAG1; RES=RESID $ CREATE; PEARSRES = RESID/SQR(STRIKES-RESID) $ CREATE; ZQ1PRED = STRIKES-RESID $ IDENTIFY; RHS = PEARSRES; pds = 15 $ IDENTIFY; RHS = PEARSRES; PDS = 6 $ ? Need to get robust standard errors CALC; NOLIST; PLOGL = LOGL; PKREG = KREG; PN = N; DEGOFFR = PN - PKREG $ MATRIX; NOLIST; PEARSON=XDOT(PEARSRES); ANB1MN = {1/DEGOFFR}*PEARSON $ ? Second, ML standard errors and t-statistics MATRIX; UNITY = INIT(PKREG,1,1) $ MATRIX; NOLIST; DIAGV = DIAG(VARB); NOLIST; PSEML = SQRT(DIAGV) | UNITY $ MATRIX; NOLIST; PTML = ISQR(DIAGV) | B $ ? Third, correct using estimated scale parameter CALC; TAU = ANB1MN $ MATRIX; LIST; TAUSQRT = SQRT(TAU); LIST; PSEMN = TAUSQRT*PSEML $ CALC; NOLIST; INVTAUSQ = 1/TAUSQRT $ MATRIX; LIST; PTMN = INVTAUSQ*PTML $ CREATE; ZQ1PRED = STRIKES-RESID $ TITLE; * RACD TABLE 7.4, col.3, p.247: Zeger-Qaqish 2 lag $ SAMPLE; 3-108 $ POISSON; LHS=STRIKES; RHS=XSTATIC,LNYSLAG1,LNYSLAG2; RES=RESID $ CREATE; PEARSRES = RESID/SQR(STRIKES-RESID) $ IDENTIFY; RHS = PEARSRES; pds = 6 $ TITLE; * RACD TABLE 7.4, col.4, p.247: Zeger-Qaqish 3 lag $ SAMPLE; 4-108 $7 POISSON; LHS=STRIKES; RHS=XSTATIC,LNYSLAG1,LNYSLAG2,LNYSLAG3; RES=RESID $ CREATE; PEARSRES = RESID/SQR(STRIKES-RESID) $ IDENTIFY; RHS = PEARSRES; pds = 6 $ ? --------- Try to estimate c in model with 1 lag TITLE; * RACD TABLE 7.4, col.5, p.247: Zeger-Qaqish 1 lag estimate c$ SAMPLE; 1-108 $ CREATE; IF (Y = 0) YSTAR1 = 1.0; (else) YSTAR1 = Y $ CREATE; LNYSTAR1 = LOG(YSTAR1) $ CREATE; IF (Y = 0) DYSTAR1 = 1.0; (else) DYSTAR1 = 0 $ CREATE; LNY1LAG1 = LNYSTAR1[-1]; DYSTLAG1 = DYSTAR1[-1] $ SAMPLE; 2-108 $ POISSON; LHS=STRIKES; RHS=XSTATIC,LNY1LAG1,DYSTLAG1; RES=RESID $ CREATE; PEARSRES = RESID/SQR(STRIKES-RESID) $ IDENTIFY; RHS = PEARSRES; PDS = 6 $ ? Need to get standard error of c using delta method TITLE; * RACD TABLE 7.4, text p.247: standard error of ZQ1 lag estimate c$ CALC; LIST; C = EXP(B(4)/B(3)); DCDB3 = C*B(4); DCDB4 = C*B(4); VARC = DCDB3*DCDB3*VARB(3,3) + DCDB4*DCDB4*VARB(4,4) + 2*DCDB3*DCDB4*VARB(3,4); SEC = SQR(VARC*213.75/102) $ ?MATRIX; DCDB = 0,0,DCDB3,DCDB4 $ ?MATRIX; LIST; DCDB $ ?MATRIX; LIST; VARC = DCDB'VARB*DCDB $ ? --------- Do NLS on Brannas model TITLE; * RACD TABLE 7.4, col.6, p.247: Brannas model NLS 1 lag $ SAMPLE; 2-108 $ NLSQ; LHS = STRIKES; LABELS = BONE,BOUTPUT,BYLAG1; START = 1.0,2.0,0.5; FCN = EXP(BONE + BOUTPUT*OUTPUT) + BYLAG1*YLAG1; pds=0 $ ? In practice need to get robust standard errors ? This will be more difficult due to nonstandard form of mean function ? This is not done here as not reported in the book ? The standard errors above assume homoskedastic error ? Get predicted mean NAMELIST; XBRANNAS = ONE,OUTPUT,YLAG1 $ CREATE; PPRED = EXP(B(1)+B(2)*OUTPUT) + B(3)*YLAG1 $ CREATE; BRANPRED = PPRED $ DSTAT; RHS = STRIKES,ZQ0PRED,ZQ1PRED,BRANPRED $ PLOT; RHS = STRIKES,ZQ0PRED,ZQ1PRED,BRANPRED $ PLOT; RHS = ZQ1PRED,BRANPRED $ STOP $