install.packages("dlm");library(dlm) install.packages("numDeriv");library(numDeriv) install.packages("tseries");library(tseries) install.packages("forecast");library(forecast) ################################## # AR(1) process ################################## # simulating the time series set.seed(4321) yt = arima.sim(n=250,list(ar=0.75,ma=0),sd=0.5) model = Arima(yt,order=c(1,0,0),include.mean=FALSE) model # set up DLM dlm0 = function(parm){ return(dlm(FF=1,V=0,GG=parm[1],W=parm[2]^2,m0=0,C0=solve(1-parm[1]^2)*parm[2]^2)) } # estimate DLM fit = dlmMLE(y=yt,parm=c(0.5,1.0),build=dlm0,hessian=T) # get estimates coef = fit$par var = solve(fit$hessian) # print results coef; sqrt(diag(var)) # get estimated variance coef[2]^2 # print DLM dlm0(fit$par) # Forecast next 5 observations using dlmForecast mod = dlm0(fit$par) modf = dlmFilter(yt,mod) fore = dlmForecast(modf,nAhead=5,method="plain") fore$f ################################## # AR(1) with Intercept ################################## set.seed(1234) n.obs = 250 yt = 2 + arima.sim(n=n.obs,list(ar=.75,ma=0),sd=.5) # set parameter restrictions (only variances here) parm_rest = function(parm){ return( c(parm[1],exp(parm[2])) ) } # set up DLM dlm1 = function(parm){ parm = parm_rest(parm) dlm = dlmModPoly(1,dV=1e-7,dW=c(0)) + dlmModARMA(ar=parm[1], ma=NULL, sigma2=parm[2]) # set initial state distribution dlm$C0[2,2] <- solve(1-parm[1]^2)*parm[2] return(dlm) } # Adding Deterministic Terms with dlmModPoly dlmModPoly(1,dV=1e-7,dW=0) # estimate parameters fit1 = dlmMLE(y=yt,parm=c(0,0),build=dlm1,hessian=T) # get parameter estimates of AR(1) part coef = parm_rest(fit1$par) # get standard errors using delta method jac = jacobian(func=parm_rest,x=fit1$par) var = jac%*%solve(fit1$hessian)%*%t(jac) # print results coef; sqrt(diag(var)) # get parameter estimates (intercept) # these are the last filtered values mod1 = dlm1(fit1$par) mod1filt = dlmFilter(yt,mod1) # get parameters coef = mod1filt$m[n.obs+1] covar = dlmSvd2var(mod1filt$U.C[[n.obs+1]], mod1filt$D.C[n.obs+1,]) coef.se = sqrt(covar[1,1]) coef; coef.se dlm1(fit1$par) # Time-varying parameter CAPM data = read.csv("sp500-ibm-monthly.csv",header=TRUE) ind = seq(1,nrow(data),length=5) date = data[ind,1] pdf(file="sp500-ibm.pdf",width=12,height=10) par(mfrow=c(2,2)) plot(data[,2],xlab="Week",ylab="Price",main="S&P500",axes=FALSE,type="l") axis(2);axis(1,at=ind,lab=date) plot(data[,3],xlab="Week",ylab="Price",main="IBM",axes=FALSE,type="l") axis(2);axis(1,at=ind,lab=date) plot(diff(log(data[,2])),xlab="Week",ylab="Log-return",main="S&P500",axes=FALSE,type="l") axis(2);axis(1,at=ind,lab=date) plot(diff(log(data[,3])),xlab="Week",ylab="Log-return",main="IBM",axes=FALSE,type="l") axis(2);axis(1,at=ind,lab=date) dev.off() r = diff(log(data[,3])) rM = diff(log(data[,2])) capm = lm(r~rM) summary(capm) pdf(file="capm.pdf",width=12,height=10) par(mfrow=c(1,1)) rang = range(rM,r) plot(rM,r,xlab="Market return (S&P500)",ylab="IBM",xlim=rang,ylim=rang) abline(capm$coef,col=2,lwd=3) title(paste("IBM = ",round(capm$coef[1],4)," + ",round(capm$coef[2],4),"*SP500",sep="")) dev.off() # set up DLM dlm2 = function(parm,x.mat){ parm = exp(parm) return( dlmModReg(X=x.mat, dV=parm[1], dW=c(parm[2],parm[3])) ) } # estimate parameters fit2 = dlmMLE(y=r,parm=c(1,1,1),x.mat=rM,build=dlm2,hessian=T) # get estimates se = sqrt(exp(fit2$par)) se # get parameter estimates over time # these are the smoothed state values mod2 = dlm2(fit2$par,rM) mod2f = dlmFilter(r,mod2) mod2s = dlmSmooth(mod2f) # plot filtered and smoothed states ind = seq(1,length(r),length=5) date = data[ind,1] pdf(file="capm-filtering-and-smoothing-intercept.pdf",width=8,height=6) plot(mod2f$m[,1],axes=FALSE,xlab="Week",ylab=expression(alpha[t]),type="l",main="") axis(2);axis(1,at=ind,lab=date) lines(mod2s$s[,1],col=2) abline(h=capm$coef[1],col=3) abline(h=0,lty=2) dev.off() pdf(file="capm-filtering-and-smoothing-slope.pdf",width=8,height=6) plot(mod2f$m[,2],axes=FALSE,xlab="Week",ylab=expression(beta[t]),type="l",main="") axis(2);axis(1,at=ind,lab=date) lines(mod2s$s[,2],col=2) abline(h=capm$coef[2],col=3) abline(h=1,lty=2) dev.off() # Quarterly UK gas consumption from 1960 to 1986, available in R as {\tt UKgas}. y = log(UKgas) pdf(file="UKgas.pdf",width=8,height=6) par(mfrow=c(1,1)) plot(log(UKgas)) dev.off() # Free-form seasonal factor DLM dlm3 = dlmModPoly() + dlmModSeas(4) buildFun = function(x) { diag(W(dlm3))[2:3] = exp(x[1:2]) V(dlm3) = exp(x[3]) return(dlm3) } fit3 = dlmMLE(y,parm=c(0.1,0.1,0.1),build=buildFun) dlm3 = buildFun(fit3$par) ySmooth = dlmSmooth(y, mod = dlm3) x = cbind(y, dropFirst(ySmooth$s[, c(1, 3)])) colnames(x) = c("Gas", "Trend", "Seasonal") pdf(file="UKgas-freeform.pdf",width=8,height=6) par(mfrow=c(1,1)) plot(x, type = "o", main = "UK Gas Consumption") dev.off() # Fourier form representation of seasonality dlm4 = dlmModPoly() + dlmModTrig(4) buildFun = function(x) { diag(W(dlm4))[2:3] = exp(x[1:2]) V(dlm4) = exp(x[3]) return(dlm4) } fit4 = dlmMLE(y,parm=c(0.1,0.1,0.1),build=buildFun) dlm4 = buildFun(fit4$par) ySmooth = dlmSmooth(y, mod = dlm4) x = cbind(y, dropFirst(ySmooth$s[, c(1, 3)])) colnames(x) = c("Gas", "Trend", "Seasonal") pdf(file="UKgas-harmonic.pdf",width=8,height=6) par(mfrow=c(1,1)) plot(x, type = "o", main = "UK Gas Consumption") dev.off()