# ClarkUCModel.ssc Clark's unobserved components model # # created: January 8, 2005 ## ## Clark model. Ref Clark (1987), JPE ## ## create state space using parameters based on estimates in MNZ (2002) ## delta = 0.8119 phi1 = 1.5303 phi2 = -0.6097 sigma.v = 0.6893 sigma.w = 0.6199 bigV = diag(c(sigma.v^2,sigma.w^2)) Omega = matrix(0,4,4) Omega[1:2,1:2] = bigV a1 = matrix(0,3,1) ## solve for initial variance of stationary part ## bigF = matrix(c(phi1,1,phi2,0),2,2) vecV = c(sigma.w^2,0,0,0) vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV P.ar2 = matrix(vecP,2,2) Sigma = matrix(0,4,3) Sigma[1,1] = -1 Sigma[2:3,2:3] = P.ar2 ## create state space list ## ssf.clark = list(mDelta=c(delta,0,0,0), mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)), mOmega=Omega, mSigma = Sigma) ## print out state space ## ssf.clark ## create state space object ## ssf.clark = CheckSsf(ssf.clark) ## simulate from model ## set.seed(569) clark.sim = SsfSim(ssf.clark,n=250) par(mfrow=c(2,1)) tsplot(clark.sim[,c(1,4)],main="Simulated response and trend") tsplot(clark.sim[,2],main="Simulated cycle") ## ## estimate Clark model for U.S. postwar quarterly real GDP ## ref: Morley, Nelson and Zivot (2002) ## ## create state space model for estimation ## Clark.mod = function(parm) { delta = parm[1] phi1 = parm[2] phi2 = parm[3] sigma2.v = exp(parm[4]) sigma2.w = exp(parm[5]) bigV = diag(c(sigma2.v,sigma2.w)) Omega = matrix(0,4,4) Omega[1:2,1:2] = bigV a1 = matrix(0,3,1) # solve for initial variance of stationary part bigF = matrix(c(phi1,1,phi2,0),2,2) vecV = c(sigma2.w,0,0,0) vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV P.ar2 = matrix(vecP,2,2) Sigma = matrix(0,4,3) Sigma[1,1] = -1 Sigma[2:3,2:3] = P.ar2 # create state space list ssf.mod = list(mDelta=c(delta,0,0,0), mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)), mOmega=Omega, mSigma = Sigma) CheckSsf(ssf.mod) } ## estimate by mle ## Clark.start=c(0.81,1.53,-0.61,-0.74,-0.96) names(Clark.start) = c("mu","phi.1","phi.2", "ln.sigma2.v","ln.sigma2.w") Clark.mle = SsfFit(Clark.start,lny.ts,"Clark.mod") ## create state space from fitted parameters ## ssf.Clark = Clark.mod(Clark.mle$parameters) ## estimated parameters and asymptotic standard errors ## summary(Clark.mle) ## estimated standard deviations and delta method standard errors ## Clark.sd = sqrt(exp(coef(Clark.mle)[4:5])) names(Clark.sd) = c("sigma.v","sigma.w") Clark.sd ## compute filtered trend and cycle ## use plot method to show filtered estimates ## filteredEst.Clark = SsfMomentEst(lny.ts,ssf.Clark,task="STFIL") names(filteredEst.Clark) colIds(filteredEst.Clark$state.moment) = c("trend(t)","cycle(t)","cycle(t-1)") plot(filteredEst.Clark, strip.text=c("trend(t)","cycle(t)","cycle(t-1)","y(t)"), main="Filtered estimates") ## create timeSeries and plot ## trend.filter = timeSeries(data=filteredEst.Clark$state.moment[,1], positions=filteredEst.Clark$positions) cycle.filter = timeSeries(data=filteredEst.Clark$state.moment[,2], positions=filteredEst.Clark$positions) par(mfrow=c(2,1)) plot(trend.filter,lny.ts, main="Log Real GDP and Filtered Trend from Clark Model", reference.grid=F) plot(cycle.filter,main="Filtered Cycle from Clark Model", reference.grid=F) ## compute smoothed trend and cycle and plot ## smoothedEst.Clark = SsfMomentEst(lny.ts,ssf.Clark,task="STSMO") trend.smooth = timeSeries(data=smoothedEst.Clark$state.moment[,1], positions=smoothedEst.Clark$positions) cycle.smooth = timeSeries(data=smoothedEst.Clark$state.moment[,2], positions=smoothedEst.Clark$positions) par(mfrow=c(2,1)) plot(trend.smooth,lny.ts, main="Log Real GDP and Smoothed Trend from Clark Model") plot(cycle.smooth,main="Smoothed Cycle from Clark Model") ## plot smoothed cycle with 95% error bands ## trend.sd = sqrt(timeSeries(data=smoothedEst.Clark$state.variance[,1], positions=smoothedEst.Clark$positions)) cycle.sd= sqrt(timeSeries(data=smoothedEst.Clark$state.variance[,2], positions=smoothedEst.Clark$positions)) upper = cycle.smooth + 2*cycle.sd lower = cycle.smooth - 2*cycle.sd par(mfrow=c(1,1)) plot(cycle.smooth,upper,lower, plot.args=list(lty=c(1,2,2)), reference.grid=F) ## Do diagnostics on standardized innovations from Kalman filter ## kf.Clark = KalmanFil(lny.ts,ssf.Clark) plot(kf.Clark) summaryStats(kf.Clark$std.innov) autocorrTest(kf.Clark$std.innov) ## KPSS specification test ## stationaryTest(lny.ts,trend="ct") ## Nyblom specification test ## n = nrow(lny.ts) ssf.Clark$mOmega[1,1] = 0 ssf.Clark$mSigma[1,1] = 0 ssf.Clark$mSigma[4,1] = smoothedEst.Clark$state.moment[1,1]- Clark.mle$parameters["mu"] kf.Clark = KalmanFil(lny.ts,ssf.Clark) test.stat = sum(cumsum(kf.Clark$innov)^2)/n^2 test.stat