## ## Nonlinear Chapter ## ## Section 2. BDS Test for Nonlinearity set.seed(10) size.mat = matrix(0, 1000, 4) for (i in 1:1000) { if (i %% 100 == 0) cat("i =", i, "\n") test.dat = rt(500, df=8) size.mat[i,] = BDSTest(test.dat, m=5, eps=1)$stat[,1] } size.p = seq(0.05, 0.95, by=0.05) size.q = qnorm(size.p) size.bds = apply(size.mat, 2, function(x) colMeans(outer(x, size.q, FUN="<="))) par(fty="s") matplot(matrix(size.p, nrow=length(size.p), ncol=4), size.bds, type="l", xlab="Nominal Size", ylab="Monte Carlo Size") legend(0.6, 0.3, paste("m=",2:5,sep=""), type="l", lty=1:4) ## BDSTest(DFX.ts, m=5) DFX.garch = garch(DFX.ts~1, ~garch(1,1), trace=F) summary(DFX.garch)$coef BDSTest(residuals(DFX.garch, standard=T), m=5, eps=c(0.5, 1, 1.5)) BDSTest(log(residuals(DFX.garch, standard=T)^2), m=5, eps=c(0.5, 1, 1.5)) ## set.seed(10) sim.garch.dat = simulate(DFX.garch, sigma=F, n.start=500, n=1000, n.rep=1000) size.garch.res = matrix(0, 1000, 4) size.garch.log = matrix(0, 1000, 4) for (i in 1:1000) { tmp = garch(sim.garch.dat[,i]~1, ~garch(1,1), trace=F) if (i %% 10 == 0) cat("Simulation No.", i, "\n") tmp.res = residuals(tmp, standardized=T) size.garch.res[i,] = BDSTest(tmp.res, m=5, eps=1)$stat[,1] size.garch.log[i,] = BDSTest(log(tmp.res^2), m=5, eps=1)$stat[,1] } size.p = seq(0.05, 0.95, by=0.05) size.q = qnorm(size.p) size.garch.res = apply(size.garch.res, 2, function(x) colMeans(outer(x, size.q, FUN="<="))) size.garch.log = apply(size.garch.log, 2, function(x) colMeans(outer(x, size.q, FUN="<="))) ## Section 3. Threshold Autoregressive Models ndx.ret2 = getReturns(ndx.dat[,"Close"])^2 ndx.rvol = sqrt(aggregate(ndx.ret2, FUN=sum, by="weeks", week.align=1)) colIds(ndx.rvol) = "RVOL" par(mfrow=c(2,2)) plot(ndx.rvol, reference.grid=F, main="RVOL") plot(log(ndx.rvol), reference.grid=F, main="Log RVOL") ndx.acf = acf(log(ndx.rvol)) ndx.pacf = acf(log(ndx.rvol), type="partial") nonlinearTest(log(ndx.rvol), method="threshold", p=6, d=1:6) nonlinearTest(log(ndx.rvol), method="threshold", p=2, d=1:2) ndx.test = nonlinearTest(log(ndx.rvol), method="threshold", p=2, d=1, save.RLS=T) par(mfrow=c(2,1)) plot(ndx.test$yd, ndx.test$tRatio[,1], xlab="Y_{t-1}", ylab="t-ratio of AR(1)") plot(ndx.test$yd, ndx.test$tRatio[,2], xlab="Y_{t-1}", ylab="t-ratio of AR(2)") ndx.setar = SETAR(log(ndx.rvol), c(-2.8, -2.4), p=2, d=1) summary(ndx.setar) par(mfrow=c(1,1)) plot(timeSeries(ndx.setar$regime, pos=positions(ndx.rvol)[-(1:2)]), reference.grid=F, ylab="Regime Index", plot.args=list(type="h")) ndx.pred = predict(ndx.setar, n.predict=100, CI.alpha=0.6, n.sim=10000) tsplot(cbind(ndx.pred$values, ndx.pred$CI), lty=c(1,6,6)) nonlinearTest(log(ndx.rvol), method="sup-LR", p=2, d=1, trim.pct=0.1, n.boot=1000) nonlinearTest(log(ndx.rvol), method="sup-LR", p=2, d=1, trim.pct=0.1, n.boot=1000, hetero=T) ndx.setar.r = TAR(log(ndx.rvol), p=2, d=1, trim.pct=0.1) ndx.setar.r$qhat summary(ndx.setar.r) names(ndx.setar.r$LR.q) plot(ndx.setar.r$LR.q$Threshold, ndx.setar.r$LR.q$LR, type="b", xlab="Threshold", ylab="LR stat") abline(h=ndx.setar.r$LR.q$Critical) ndx.pred.2 = predict(ndx.setar.r, n.predict=100, CI.alpha=0.6, n.sim=10000) nonlinearTest(log(ndx.rvol), method="STAR-LM", p=2, d=1:2) nonlinearTest(log(ndx.rvol), method="STAR-LM", p=2, d=1:2, hetero=T) ## Section 4. Smooth Transition Autoregressive Models ndx.lstar = STAR(log(ndx.rvol), p=2, d=1) summary(ndx.lstar) ndx.pred.3 = predict(ndx.lstar, n.predict=100, CI.alpha=0.6, n.sim=10000) tsplot(cbind(ndx.pred.3$values, ndx.pred.3$CI), lty=c(1,6,6)) ESTAR.res = function(theta, g.scale, x, y, q) { k = ncol(x) G = 1 - exp(-exp(theta[1])/g.scale * (q - theta[2])^2) X = cbind(x * (1 - G), x * G) m = crossprod(t(backsolve(chol(crossprod(X)), diag(2 * k)))) beta = m %*% t(X) %*% y y - X %*% beta } ndx.LHS = log(ndx.rvol)[3:length(ndx.rvol)]@data ndx.RHS = cbind(1, tslag(log(ndx.rvol), 1:2, trim=T)@data) ndx.estar = nlregb(length(ndx.rvol)-2, start=c(0,mean(ndx.RHS[,2])), residuals=ESTAR.res, lower=c(-Inf, min(ndx.RHS[,2])), upper=c( Inf, max(ndx.RHS[,2])), g.scale=var(ndx.RHS[,2]), x=ndx.RHS, y=ndx.LHS, q=ndx.RHS[,2]) ndx.estar$parameters exp(ndx.estar$parameters[1])/var(ndx.RHS[,2]) ## Section 5. Markov Switching State Space Models mchain.p0(matrix(c(0.47, 0.05, 0.53, 0.95), 2, 2)) GetSsfMSAR = function(parm) { mDelta = mDelta.other = rep(0, 3) mDelta[1] = parm[1] mDelta.other[1] = parm[2] # mPhi = mPhi.other = matrix(0, 3, 2) mPhi[1,] = c(parm[3], parm[4]) mPhi.other[1,] = c(parm[5], parm[6]) mPhi[2:3,1] = mPhi.other[2:3,1] = 1 # mOmega = mOmega.other = matrix(0, 3, 3) mOmega[1,1] = parm[7] mOmega.other[1,1] = parm[8] # mSigma = matrix(0, 3, 2) mSigma[1:2, 1:2] = diag(1e+6, 2) # mTrans = matrix(0, 2, 2) mTrans[1,1] = parm[9] mTrans[1,2] = 1 - mTrans[1,1] mTrans[2,2] = parm[10] mTrans[2,1] = 1 - mTrans[2,2] # list(mDelta=mDelta, mDelta.other=mDelta.other, mPhi=mPhi, mPhi.other=mPhi.other, mOmega=mOmega, mOmega.other=mOmega.other, mSigma=mSigma, mTrans=mTrans) } GetSsfMSAR2 = function(parm) { parm = as.vector(parm) # mDelta = mDelta.other = rep(0, 3) mDelta[1] = parm[1] mDelta.other[1] = parm[1] + exp(parm[2]) # AR11 = 2/(1+exp(-parm[3])) - 1 AR12 = 2/(1+exp(-(parm[3]+exp(parm[4])))) - 1 AR21 = 2/(1+exp(-parm[5])) - 1 AR22 = 2/(1+exp(-(parm[5]+exp(parm[6])))) - 1 # mPhi = mPhi.other = matrix(0, 3, 2) mPhi[1,] = c(AR11+AR12, -AR11*AR12) mPhi.other[1,] = c(AR21+AR22, -AR21*AR22) mPhi[2:3,1] = mPhi.other[2:3,1] = 1 # mOmega = matrix(0, 3, 3) mOmega[1,1] = exp(parm[7]) # mSigma = matrix(0, 3, 2) mSigma[1:2, 1:2] = diag(1e+6, 2) # mTrans = matrix(0, 2, 2) mTrans[1,2] = 1/(1+exp(-parm[8])) mTrans[1,1] = 1 - mTrans[1,2] mTrans[2,1] = 1/(1+exp(-parm[9])) mTrans[2,2] = 1 - mTrans[2,1] # ssf = list(mDelta=mDelta, mDelta.other=mDelta.other, mPhi=mPhi, mPhi.other=mPhi.other, mOmega=mOmega, mTrans=mTrans, mSigma=mSigma) CheckSsf(ssf) } ndx.start = c(-2, -0.7, -0.7, 0.7, -0.7, 0.7, -2, -2, -3) names(ndx.start) = c("mu1", "mu2", "phi11", "phi12", "phi21", "phi22", "sigma", "p", "q") ndx.msar = SsfFitMS(ndx.start, log(ndx.rvol), GetSsfMSAR2, l.start=11) class(ndx.msar) summary(ndx.msar) ndx.ssf = GetSsfMSAR2(ndx.msar$parameters) cbind(ndx.ssf$mDelta, ndx.ssf$mDelta.other) ndx.ssf$mPhi ndx.ssf$mPhi.other ndx.ssf$mOmega ndx.ssf$mTrans ndx.f = SsfLoglikeMS(log(ndx.rvol), ndx.ssf, save.rgm=T) par(mfrow=c(2,1)) plot(timeSeries(ndx.f$regimes[,1], pos=positions(ndx.rvol)), reference.grid=F, main="Filtered Low Vol Regime Prob") plot(timeSeries(ndx.f$regimes[,3], pos=positions(ndx.rvol)), reference.grid=F, main="Smoothed Low Vol Regime Prob") ## Section 6. An Extended Example GetSsfCoinIndex = function(parm) { parm = as.vector(parm) mDelta = mDelta.other = rep(0, 11) mDelta[1] = parm[1] mDelta.other[1] = parm[1] + exp(parm[2]) # AR.C1 = 2/(1+exp(-parm[3])) - 1 AR.C2 = 2/(1+exp(-(parm[3]+exp(parm[4])))) - 1 # AR.e1 = 2/(1+exp(-parm[5])) - 1 AR.e2 = 2/(1+exp(-parm[6])) - 1 AR.e3 = 2/(1+exp(-parm[7])) - 1 AR.e4 = 2/(1+exp(-parm[8])) - 1 # mPhi = matrix(0, 11, 7) mPhi[1,1:2] = c(AR.C1+AR.C2, -AR.C1*AR.C2) mPhi[2,1] = 1 mPhi[3,3] = AR.e1 mPhi[4,4] = AR.e2 mPhi[5,5] = AR.e3 mPhi[6,6] = AR.e4 mPhi[7,1] = mPhi[7,7] = 1 # mPhi[8:10,1] = parm[9:11] mPhi[11,1:2] = parm[12:13] mPhi[8,3] = mPhi[9,4] = mPhi[10,5] = mPhi[11,6] = 1 # mOmega = matrix(0, 11, 11) mOmega[1,1] = 1 mOmega[3,3] = exp(parm[14]) mOmega[4,4] = exp(parm[15]) mOmega[5,5] = exp(parm[16]) mOmega[6,6] = exp(parm[17]) # mTrans = matrix(0, 2, 2) mTrans[1,2] = 1/(1+exp(-parm[18])) mTrans[1,1] = 1-mTrans[1,2] mTrans[2,1] = 1/(1+exp(-parm[19])) mTrans[2,2] = 1-mTrans[2,1] # mSigma = matrix(0, 8, 7) mSigma[1:7, 1:7] = diag(1e+6, 7) ans = list(mDelta=mDelta, mDelta.other=mDelta.other, mSigma=mSigma, mOmega=mOmega, mPhi=mPhi, mTrans=mTrans) CheckSsf(ans) } DOC.dat = getReturns(DOC.ts[,1:4], percentage=T) DOC.dat@data = t(t(DOC.dat@data) - colMeans(DOC.dat@data)) DOC.dat@data = t(t(DOC.dat@data) / colStdevs(DOC.dat@data)) DOC.start = c(-1.5, 0.6, 0.3, 0.1, .1, .1, .1, .1, 0.3, 0.3, 0.3, 0.3, 0.1, -.5, -.5, -.5, -.5, -1.5, -3) names(DOC.start) = c("mu1", "mu2", "phi1", "phi2", "psi1", "psi2", "psi3", "psi4", "L1", "L2", "L3", "L41", "L42", "s1", "s2", "s3", "s4", "p", "q") DOC.fit = SsfFitMS(DOC.start, DOC.dat, GetSsfCoinIndex, l.start=13, trace=T) summary(DOC.fit) DOC.ssf = GetSsfCoinIndex(DOC.fit$parameters) c(DOC.ssf$mDelta[1], DOC.ssf$mDelta.other[1]) print(DOC.ssf$mPhi, digits=3) DOC.f = SsfLoglikeMS(DOC.dat, DOC.ssf, save.rgm=T, l.start=13) names(DOC.f) DOC.dates = positions(DOC.dat)[-(1:12)] filt.p = timeSeries(DOC.f$regimes[,1], pos=DOC.dates) smoo.p = timeSeries(DOC.f$regimes[,3], pos=DOC.dates) par(mfrow=c(2,1)) plot(filt.p, reference.grid=F, main="Filtered Recession Probability") plot(smoo.p, reference.grid=F, main="Smoothed Recession Probability") DOC.index = lm(DOC.ts@data[,5]~I(1:433))$residuals[-(1:13)] filt.ci = rowSums(DOC.f$state[,c(7,14)]*DOC.f$regime[,1:2]) filt.ci = timeSeries(filt.ci, pos=DOC.dates) plot(filt.ci, reference.grid=F, main="Filtered MS Coincident Index") doc.ci = timeSeries(DOC.index, pos=DOC.dates) plot(doc.ci, reference.grid=F, main="DOC Coincident Index") # Use Andrew's shadePlot to add business cycle dates to plot NBERShade = list(x=data.frame(timeDate(c("04/01/60","10/01/69","11/01/73","01/01/80","07/01/81","07/01/90")), timeDate(c("02/01/61","11/01/70","03/01/75","07/01/80","11/01/82","03/01/91"))), col=12) par(mfrow=c(2,1)) shadePlot(filt.p, reference.grid=F, main="Filtered Recession Probability", shade.regions=list(NBERShade)) shadePlot(smoo.p, reference.grid=F, main="Smoothed Recession Probability", shade.regions=list(NBERShade))