## StAR Model#####
## Random number seed
set.seed(4093)
## Creating trend variable.
t <- seq(1,50,1)
# Generating data on y.
y <- 0.004 + 0.0045*t - 0.09*t^2 + 50*rt(50,df=5)
T <- length(y)
# The trend matrix
Trend <- cbind(1,poly(t,2,raw=TRUE))
# Estimating the model
star <- StAR(y,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.
oldpar <- par(mfcol=c(3,1))
matplot(dates[2:T],cbind(y[2:T],star$fitted,star$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(star$res,main="",xlab="") ## Histogram of y
matplot(dates[2:T],cbind(star$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