## StVAR Model#####
## Random number seed
set.seed(7504)
## Creating trend variable.
t <- seq(1,50,1)
y <- x <- vector(length=50)
y[1] <- x[1] <- 0
# Generating data on y and x.
ut <- rt(50,df=5)
for(i in 2:50)
{
y[i] <- 0.004 + 0.0045*t[i] - 0.9*y[i-1] + .2*x[i-1] + ut[i]
x[i] <- 0.05 - 0.005*t[i] + 0.8*x[i-1] + 1*y[i-1] + ut[i]
}
# The trend matrix
Trend <- cbind(1,t)
# Estimating the model
stvar <- StVAR(cbind(y,x),lag=1,Trend=Trend,v=6,maxiter=2000,hes=TRUE)
# 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:length(y)],cbind(y[2:length(y)],stvar$fit[,1],stvar$trend[,1]),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(stvar$res[,1],main="",xlab="") ## Histogram of y
matplot(dates[2:length(y)],cbind(stvar$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)
## Plotting the variable x, its estimated trend and the fitted value.
oldpar <- par(mfcol=c(3,1))
matplot(dates[2:length(x)],cbind(x[2:length(x)],stvar$fit[,2],stvar$trend[,2]),
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(stvar$res[,2],main="",xlab="") ## Histogram of x
matplot(dates[2:length(x)],cbind(stvar$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