#PROGRAM 5 - SIMULATING FROM A SPECIFIED MODEL AND THEN FORECASTING
#
# K. Ensor, ensor@rice.edu
# Rice University
#

# Let's first create a time series of length nobs.
# Just as an example, I have defined daily data that
# has a weekly period.
#
# Change the model as you would like by changing the
# set of specifications below. The final list, kmodel
# is the one that feeds into the other routines.

# Plotting function used later

prettyplot<-function(nrow,ncol){
	par(mfrow=c(nrow,ncol),omi=c(.5,1,.5,1))
}

#define model, change as desired
kmodel1<-list(ar=c(.4),period=1,sigma2=4) #time unit one day
kmodel2<-list(ar=c(.7),period=5,sigma2=.5) #time unit one week
kmodel<-list(kmodel1,kmodel2)
kreg=rep(1,nobs)		#vector of 1's for the intercept
kintercept=25			#intercept for additive regression component

#number of observations and number of forecasts
nobs=200
nfore=10
k=seq((nobs-50),nobs) # sets the number of observations to plot 
						  # before the forecast
						  # to get all observations in the series
						  # change (nobs-50) to 1

#simulate an arima model with additive regression variables
x=arima.sim(nobs,model=kmodel,xreg=kreg,reg.coef=kintercept)

prettyplot(2,2)  # defined in MISC
# examine the basic plots for the simualted series
tsplot(x)
hist(x)
acf(x)
pacf(x)			

# Fit the model and forecast from the fitted model
kmodel.fit=arima.mle(x,kmodel,xreg=1)
kmodel.fore<-arima.forecast(x,n=nfore,model=kmodel.fit$model,
	xreg=rep(1,(nobs+nfore)),reg.coef=kmodel.fit$reg.coef)
kerr=1.96*kmodel.fore$std.err
kf=kmodel.fore$mean

prettyplot(1,1)
#compute the mean, variance and standard error of the mean
xm=mean(x)
xmvar=var(x)/nobs
xmv=rep(xm,length(kf))
xmsdr=rep((1.96*(xmvar^.5)),length(kf))
xmvar=var(x)
xsdr=rep((1.96*(xmvar^.5)),length(kf))
						
#plot the model based and iid prediction intervals 
tsplot(c(x[k],xmv),
		c(x[k],(xmv-xsdr)),
		c(x[k],(xmv+xsdr)),   #standard CI
 		c(x[k],kf),
		c(x[k],(kf-kerr)),
		c(x[k],(kf+kerr)),
		x[k],
		rep(kintercept,(length(k)+nfore)),
		col=c(6,6,6,3,3,3,1,1),
		lty=c(1,1,1,1,1,1,1,1)) 		 # ts pred int
title(main="Simulated Series + Model Based Forecast and 95% PB \n+ Sample Mean and 95% PI")

#plot the model based prediction intervals and the iid CI
tsplot(c(x[k],xmv),
		c(x[k],(xmv-xmsdr)),
		c(x[k],(xmv+xmsdr)),   #standard CI
 		c(x[k],kf),
		c(x[k],(kf-kerr)),
		c(x[k],(kf+kerr)),
		x[k],
		rep(kintercept,(length(k)+nfore)),
	   col=c(6,6,6,3,3,3,1,1),
	   lty=c(1,1,1,1,1,1,1,1)) 		 # ts pred int
title(main="Simulated Series + Model Based Forecast and 95% PB \n+ Sample Mean and 95% CI")