SPECIFIC R COMMANDS FOR TIME SERIES ANALYSIS

postscript(file="plot.ps",horizontal=F)    ## Plots directed to file "plot.ps"
par(mfrow=c(3,1))

## Basic Time Series Plot, Sample ACF and PACF

yt <- scan("/u/r/e/reinsel/stat349/grind.dat")    ## Read in time series data
ny <- length(yt)
library(ts)                                ## Loads the library for time series
plot.ts(yt, type="b", pch="*", main="Grinding Wheel Profile Data")
yt.acr <- acf(yt, lag.max=25, plot = F)             ## Create ACF for series yt
plot.acf(yt.acr, main="ACF of Original Series yt")
acf(yt,xlim=c(2,26),ylim=c(-0.75,0.75),main="ACF without lag-0 term")
cat("\n ACF for original series yt",fill=T)
print(yt.acr)
yt.pacf <- acf(yt, lag.max=25, type="partial", plot = F)  ## Create PACF for yt
plot.acf(yt.pacf, main="PACF of Original Series yt")
cat("\n PACF for original series yt",fill=T)
print(yt.pacf)

wt <- diff(yt)                                ## Create 1st differences of yt

## Specifying ARIMA Models, Fitting, Model Checking, and Forecasting

## General form of multiplicative seasonal ARIMA model statement
#arima(yt, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S), method="ML")
## Estimate ARIMA model with additional regressor terms x1 and x2 included
#yfit1 <- arima(yt, order=c(3,0,0), method="ML", xreg=cbind(x1,x2)) 

par(mfrow=c(1,1))
yfit1 <- arima(yt, order=c(3,0,0), method="ML")   ## Estimate ARIMA model by ML
cat("\n ARIMA(3,0,0) model for original series yt",fill=T)
print(yfit1)
diag1 <- tsdiag(yfit1)                  ## Create Diagnostics for fitted model
res.acr <- acf(resid(yfit1), lag.max=25, plot = F)  ## Create ACF for residuals
cat("\n ACF for residuals from ARIMA(3,0,0) model",fill=T)
print(res.acr)
cat("\n Box-Pierce statistics of residuals from ARIMA(3,0,0) model",fill=T)
Box.test(resid(yfit1), lag = 12, type = "Ljung-Box")  ##Obtain Ljung-Box stat.
forecast1 <- predict.Arima(yfit1, n.ahead=10, newreg=10)  ## Obtain forecasts for model        
cat("\n Forecasts of yt for 10 time periods ahead based on ARIMA(3,0,0) model",fill=T)
print(forecast1)

q()                                               ## Ends R session