### ARIMA simulations, estimation, forecasting

### Note that in Splus the ARIMA model is defined as follows.

### x_t = phi_1 x_t-1 + ... + phi_p x_t-p + w_t - theta_1 w_t-1 - ... - theta_q w_t-q

### This is different from the book in the MA coefficients.  In Splus the coefficents 
### have negative signs and in the book they have positive signs.

######################################################################################
# Simulation

# MA(1)

par(mfrow=c(2,3))

x <- arima.sim(100,model=list(ma=-.9))	# for the book theta = +.9
x1 <- x
ts.plot(x,main="MA(1), theta=.9")
acf(x, 20)
acf(x, 20, "partial")

x <- arima.sim(100,model=list(ma=.9))		# for the book theta = -.9
x2 <- x
ts.plot(x,main="MA(1), theta=-.9")
acf(x, 20)
acf(x, 20, "partial")

# MA(2)

par(mfrow=c(2,3))

x <- arima.sim(400,model=list(ma=c(-.6,-.5))) 	# for the book theta_1 = .6, theta_2 = .5
x3 <- x
tsplot(x, main="MA(2),theta_1=.6, theta_2=.5")
acf(x, 20)
acf(x, 20, "partial")

x <- arima.sim(400,model=list(ma=c(.6,.5)))		# for the book theta_1 = -.6, theta_2 = -.5
x4 <- x
tsplot(x, main="MA(2),theta_1=-.6, theta_2=-.5")
acf(x, 20)
acf(x, 20, "partial")

par(mfrow=c(2,3))

x <- arima.sim(400,model=list(ma=c(-.6,.5)))		# for the book theta_1 = .6, theta_2 = -.5
x5 <- x
tsplot(x, main="MA(2),theta_1=.6, theta_2=-.5")
acf(x, 20)
acf(x, 20, "partial")

x <- arima.sim(400,model=list(ma=c(.6,-.5)))		# for the book theta_1 = -.6, theta_2 = .5
x6 <- x
tsplot(x, main="MA(2),theta_1=-.6, theta_2=.5")
acf(x, 20)
acf(x, 20, "partial")

######################################################################################

# AR(1)

par(mfrow=c(2,3))

x <- arima.sim(200,model=list(ar=.5))
x7 <- x
tsplot(x, main="AR(1),phi=.5")
acf(x, 20)
acf(x, 20, "partial")

x <- arima.sim(200,model=list(ar=-.5))
x8 <- x
tsplot(x,main="AR(1),phi=-.5")
acf(x, 20)
acf(x, 20, "partial")

# AR(2)

par(mfrow=c(2,3))

x <- arima.sim(400,model=list(ar=c(.4,.2)))
x9 <- x
tsplot(x,main="AR(2),phi_1=.4,phi_2=.2")

acf(x, 20)
acf(x, 20, "partial")

# AR(2) with complex roots, see Example 2.9 page 105

x <- arima.sim(144,model=list(ar=c(1.5,-.75)))
x10 <- x
tsplot(x,main="AR(2),phi_1=1.5,phi_2=-.75, complex roots")
acf(x, 20)
acf(x, 20, "partial")

######################################################################################

# ARMA(1,1)

par(mfrow=c(2,3))

x <- arima.sim(100,model=list(ar=.5,ma=-.6))
x11 <- x
tsplot(x,main="ARMA(1,1),phi=.5,theta=-.6")
acf(x, 20)
acf(x, 20, "partial")

# ARIMA(1,1,1)

x <- arima.sim(100,model=list(ar=.5,ndiff=1,ma=-.6))	# nonstationary
x12 <- x
tsplot(x,main="AR1MA(1,1,1),phi=.5,diff=1,theta=-.6")
acf(x, 20)
acf(x, 20, "partial")

######################################################################################
# Estimation

# Fit an AR(p) model to x7 		Ans: p = 1

x7.ar <- ar(x7, order.max=5)		# order.max determines the maximum number of lags p
										# to try 
										# the default is method="yule-walker"
										
x7.ar		# output

par(mfrow=c(2,3))

x <- 0:5								# plot the aic for models p = 0,1,...,p
plot(x,x7.ar$aic,type="l")
ts.plot(x7.ar$resid)					# plot residuals, should be while noise if model fits
acf(x7.ar$resid, 20)					# ACF should have no significant spikes
acf(x7.ar$resid, 20, "partial")	# PACF should have no significant spikes
hist(x7.ar$resid)						# the residuals should be normal

# Fit an AR(p) model to x9 		Ans: p = 2

x9.ar <- ar(x9, order.max=5, method="burg")	
										
x9.ar	

par(mfrow=c(2,3))

x <- 0:5							
plot(x,x9.ar$aic,type="l")
ts.plot(x9.ar$resid)				
acf(x9.ar$resid, 20)				
acf(x9.ar$resid, 20, "partial")
hist(x9.ar$resid)					

# Fit an MA(q) model to x2		Ans: q = 1, theta = -.9

x2.arma <- arima.mle(x2,list(order=c(0,0,1)))

x2.arma

# Fit an ARMA(1,1) to x11

x11.arma <- arima.mle(x11,list(order=c(1,0,1)))

x11.arma

# Fit an ARIMA(1,1,1) to x12

x12.arma <- arima.mle(x12,list(order=c(1,1,1)))

x12.arma