### Statistics 6871 - Time Series

### Lab 6

### A few more examples of fitting ARIMA models.

### All of the data sets are from Pankratz, Forecasting With Univariate Box-Jenkins Models

### Again note that Splus is case sensitive!

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

# QUARTERLY CHANGE IN BUSINESS INVENTORIES (BILLIONS) 1955-1969

# In this example we analyze the quarterly change in business inventories, stated at annual

# rates in billions of dollars. We examime 60 observations covering the period from the first

# quarter 1955 through the fourth quarter of 1969.

# That this data has already been seasonally adjusted.

# Load the data

It <- rts(scan("I:\\Courswrk\\Stat\\esuess\\Stat6871\\Data_p\\Inventory.txt"))

# Subtract out the mean.

It <- It - mean(It)

# Plot the data. We might say the data looks aproximately stationary, remember we are only

# looking at 60 data points.

par(mfcol=c(3,1))

ts.plot(It)

# Plot the ACF and PACF with the idea of identifying a possible ARMA model for the data.

# With so few data points at rule of n/4 is used when determining the number of lags to plot

# in the ACF. So for this set of data we use lag = 60/4 =15. Recall that the AR model is

# suggested by the PACF and the MA model is suggested by the ACF. For identification purposes

# we want to find a model with the least number of coeficients to estimate.

acf(It, lag = 15)

acf(It, lag = 15, type = "partial")

# From the PACF we think that an AR model might be a good chose to model the It series.

# Fit an AR(1) and record the paramter estimate and wither it is statistically significant.

model <- list(order = c(1,0,0))

It.arima100 <- arima.mle(It, model)

It.arima100$model

# Diagnostics.

par(mfrow = c(1,1))

It.arima100.diag <- arima.diag(It.arima100)

# Forcasts. Suppose we are interested in forecasting the next 4 quarters. Splus has a

# command that will produce such forecasts. The command is arima.forecast, it has three

# arguments, the time series data, the model you got from fitting the ARIMA model, and the

# number of forecasts you want to make.

It.fore <- arima.forecast(It, model = It.arima100$model, n = 4)

# We can also plot a picture with the data and forecasts with approximate confidence intervals.

tsplot(It,It.fore$mean,It.fore$mean+2*It.fore$std.err,It.fore$mean-2*It.fore$std.err)

# What do you notice about the forecasting intervals?

# Next I would like you try the same analysis in Minitab. Load the data set Inventory1.txt

# into Minitab using the Other Files, Special text command under the option File. Under the

# Stat option in Minitab there are some Time Series functions. Go through the above

# modelling and compare your answers. When you go to fit an ARIMA model go into the Graphs

# option and click on everything and then go into the Forecast option and give it a lead of 4.

# Be sure to got to the help file and read the page Time Series Overview.

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

 

# PERSONAL SAVING AS % OF DISPOSABLE INCOME 1955-1979

# Try the same kind of modelling with the following data set. In this example we will analyze

# 100 quarterly observations of the U.S. savsings rate for the years 1955-1979. The data are

# seasonally adjusted prior to publication by the U.S. Department of Commerce. We will assume

# the data is approximately stationary to begin with.

St <- rts(scan("I:\\Courswrk\\Stat\\esuess\\Stat6871\\Data_p\\SaveRate.txt"))

ts.plot(St)

acf(St)

acf(St, type="partial")

# Continue. Make sure you examine any residuals you come up with closely.

# Again try this data analysis in Minitab. Load in SavsRate1.txt.

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

# MONTHLY BITUMINOUS COAL PRODUCTION IN U.S.A. 1952-1959

Coal <- rts(scan("I:\\Courswrk\\Stat\\esuess\\Stat6871\\Data_p\\Coal.txt"))

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

# NEW PRIVATE HOUSING PERMITS ( QUARTERLY 1947-1967 )

Permits <- rts(scan("I:\\Courswrk\\Stat\\esuess\\Stat6871\\Data_p\\Permits.txt"))