## StDLM Model#####
## Random number seed
set.seed(7504)
## Creating trend variable.
t <- seq(1,100,1)
ut <- rt(100,df=5)
# Generating data on y, x and z.
y <- 0.004 + 0.0045*t - 0.09*t^2 + 0.001*t^3 + 50*ut
x <- 0.05 - 0.005*t + 0.09*t^2 - 0.001*t^3 + 40*ut
z <- 0.08 - 0.006*t + 0.08*t^2 - 0.001*t^3 + 30*ut
# The trend matrix
Trend <- cbind(1,poly(t,3,raw=TRUE))
# Estimating the model
stdlrm <- StDLM(y,cbind(x,z),lag=1,Trend=Trend,v=5,maxiter=2000)
# Generate arbitrary dates
dates <- seq(as.Date("2014/1/1"), as.Date("2016/1/1"), "weeks")
## Plotting the variable y, its estimated trend and the fitted value.
lag <- 1
oldpar <- par(mfcol=c(3,1))
matplot(dates[(lag+1):length(y)],cbind(y[(lag+1):length(y)],stdlrm$fit,stdlrm$trend),
xlab="Months",type='l',lty=c(1,2,3),lwd=c(1,1,3),col=c("black","blue","black"),ylab=" ",xaxt="n")
axis.Date(1, at = seq(as.Date("2014/1/1"), as.Date("2016/1/1"), "months"),labels=TRUE)
legend("bottomleft",legend=c("data","trend","fitted values"),lty=c(1,2,3),lwd=c(1,1,3),
col=c("black","blue","black"),cex=.85)
hist(stdlrm$res,main="",xlab="") ## Histogram of y
matplot(dates[2:length(y)],cbind(stdlrm$cvar),xlab="Months",type='l',lty=2,lwd=1,
ylab="fitted variance",xaxt="n")
axis.Date(1, at = seq(as.Date("2014/1/1"), as.Date("2016/1/1"), "months"),labels=TRUE)
par(oldpar)
Run the code above in your browser using DataLab