## stlm Model#####
## Random number seed
set.seed(7504)
## Creating trend variable.
t <- seq(1,50,1)
# Generating data on y, x and z.
x <- rt(50,df=5) ; set.seed(7505)
z <- 5*rt(50,df=5) ; set.seed(7506)
y <- 4 + 0.045*t + 2*x + 5*z + rt(50,25,df=5)
# The trend matrix
Trend <- cbind(1,poly(t,1,raw=FALSE))
colnames(Trend) <- c("const","t")
# Estimating the model
stlrm <- StLM(y,cbind(x,z),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.
oldpar <- par(mfcol=c(3,1))
matplot(dates[1:length(y)],cbind(y[1:length(y)],stlrm$fit,stlrm$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(stlrm$res,main="",xlab="") ## Histogram of y
matplot(dates[1:length(y)],cbind(stlrm$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