# MNZModel.ssc Morley, Nelson and Zivot's UC model # # Created: January 8, 2005 ## ## MNZ model: Ref. Morley, Nelson and Zivot (2002), "Why are ## Beveridge-Nelson and Unobserved Components Decompositions of GDP ## so Different?", forthcoming in Review of Economics and Statistics ## ## Note: parameters are taken from Table 4 of ## Morley, Nelson and Zivot (2002) ## delta = 0.8156 phi1 = 1.3419 phi2 = -0.7060 sigma.v = 1.2368 sigma.w = 0.7485 sigma.vw = -0.8389 bigV = matrix(c(sigma.v^2,sigma.vw,sigma.vw,sigma.w^2),2,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 ssf.mnz= 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) ssf.mnz ## simulate from model ## set.seed(569) mnz.sim = SsfSim(ssf.mnz,n=250) par(mfrow=c(2,1)) tsplot(mnz.sim[,c(1,4)],main="Simulated response and trend") tsplot(mnz.sim[,2],main="Simulated cycle") ## estimate model using US postwar quarterly real GDP ## MNZ.mod = function(parm) { delta = parm[1] phi1 = parm[2] phi2 = parm[3] sigma.v = exp(parm[4]) sigma.w = exp(parm[5]) rho.vw = parm[6] sigma.vw = sigma.v*sigma.w*rho.vw bigV = matrix(c(sigma.v^2,sigma.vw,sigma.vw,sigma.w^2),2,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 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) } ## set starting values ## MNZ.start=c(0.81,1.34,-0.70,0.21,-0.30,-0.9) names(MNZ.start) = c("delta","phi.1","phi.2", "ln.sigma.v","ln.sigma.w","rho") ## estimate model by MLE ## low.vals = c(0,0,-2,-Inf,-Inf,-0.999) up.vals = c(2,2,2,Inf,Inf,0.999) MNZ.mle = SsfFit(MNZ.start,lny.ts,"MNZ.mod", lower=low.vals,upper=up.vals) ## estimated parameters and asymptotic standard errors ## MNZ.mle$parameters sqrt(diag(MNZ.mle$vcov)) # estimated standard deviations exp(MNZ.mle$parameters[4:5]) # note: use delta method to get standard errors for standard deviations # compute filtered trend and cycle filteredEst.MNZ = SsfMomentEst(lny.ts, MNZ.mod(MNZ.mle$parameters),task="STFIL") trend.filter = timeSeries(data=filteredEst.MNZ$state.moment[,1], positions=filteredEst.MNZ$positions) cycle.filter = timeSeries(data=filteredEst.MNZ$state.moment[,2], positions=filteredEst.MNZ$positions) par(mfrow=c(2,1)) plot(trend.filter,lny.ts, main="Log Real GDP and Filtered Trend from MNZ Model",reference.grid=F) plot(cycle.filter,main="Filtered Cycle from MNZ Model",reference.grid=F) # compute smoothed trend and cycle smoothedEst.MNZ = SsfMomentEst(lny.ts, MNZ.mod(MNZ.mle$parameters),task="STSMO") trend.smooth = timeSeries(data=smoothedEst.MNZ$state.moment[,1], positions=smoothedEst.MNZ$positions) cycle.smooth = timeSeries(data=smoothedEst.MNZ$state.moment[,2], positions=smoothedEst.MNZ$positions) par(mfrow=c(2,1)) plot(trend.smooth,lny.ts, main="Log Real GDP and Smoothed Trend from MNZ Model") plot(cycle.smooth,main="Smoothed Cycle from MNZ Model") plot(cycle.smooth,reference.grid=F) # plot smoothed cycle with 95% error bands cycle.sd= sqrt(timeSeries(data=smoothedEst.MNZ$state.variance[,2], positions=smoothedEst.MNZ$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,3,3))) # plot filtered and smoothed cycle together plot(cycle.filter,cycle.smooth, plot.args=list(lty=c(1,2))) legend(0.75,-3,legend=c("filtered","smoothed"),lty=c(1,2)) # plot filtered and smoothed trend together plot(diff(trend.filter),diff(trend.smooth))