Advanced Statistics with R code 90% doen just check it and make a summary

CSC425 – Time series analysis and forecasting

Homework 4 - Solutions

Problem 1 [just for fun!] – NO CREDIT

I noticed some students had questions regarding the SARIMA notation. Write the model expression for

the following SARIMA models for the {Xt} time series.

Example: SARIMA(2,0,0)(1,0,0)[12]

The notation indicates a regular AR(2) model and a seasonal AR(1)[12] model

(1-φ1B-φ2B 2 )(1-Φ1 B

12 )Xt = at

a) SARIMA(2,1,0)(1,1,0)[4] b) SARIMA(1,0,1)(0,0,1)[12] c) SARIMA(0,1,3)(0,1,1)[4]

In general SARIMA(p,d,q)(P,D,Q)[s] where the first set (p,d,q) describes the regular

ARIMA(p,d,q) model, and the second set (P,D,Q)[s] describes the order of the seasonal

component.

a) The first set (2,1,0) identifies a first difference + AR(2) model. The second set (1,1,0) indicates the seasonal component with a seasonal differencing and a AR(4) seasonal

model. Thus

SARIMA(2,1,0)(1,1,0)[4]

(1-ϕ1B-ϕ2B 2 )(1-φ1B

4 ) (1-B)(1-B

4 )Xt = at

b) The first set (1,0,1) indicates an ARMA(1,1) model. The second set (0,0,1) is for the seasonal component, indicating MA(12) seasonal model. Thus

SARIMA(1,0,1)(0,0,1)[12]

(1-ϕ1B)Xt =(1-θ1B 12

) (1-ϑ1B)at

c) The first set (0,1,3) indicates an ARIMA(0,1,3) model, i.e. first difference + MA(3) model. The second set (0,1,1,) indicates a seasonal differencing and a MA(4) model.

Thus

SARIMA(0,1,3)(0,1,1)[4]

(1-B)(1-B 4 )Xt = (1-θ1B

4 ) (1-ϑ1B- ϑ2B

2 - ϑ3B

3 )at

Problem 2 [15pts]

Consider the quarterly earnings per share (EPS) of Costco from the last quarter of 2001 to the third

quarter of 2013 fiscal year. The data were obtained from the Nasdaq website and are saved in the

costco.csv file.

a) Analyze the distribution of the quarterly earnings. Are the data normally distributed? Provide appropriate statistics, tests and graphs to support your conclusions.

The average quarterly earnings are $0.66 per share with a standard deviation of $0.27. The

distribution of eps has a slight positive skew – as shown by the histogram and normal quantile plot

(points under the right tails diverge from the diagonal). Normality test rejects hypothesis of

normality for earnings data at 5% level. (See SAS / R output)

The UNIVARIATE Procedure

Variable: eps

Moments

N 48 Sum Weights 48

Mean 0.66270833 Sum Observations 31.81

Std Deviation 0.26970228 Variance 0.07273932

Skewness 0.95652585 Kurtosis 0.80234438

Uncorrected SS 24.4995 Corrected SS 3.41874792

Coeff Variation 40.6969796 Std Error Mean 0.03892817

Results of Jacque and Bera test on normality

Jarque & P-value for

skewness, kurtosis, Bera Jarque &

Obs eps eps statistic Bera test

1 0.95653 0.80234 8.60705 0.013521

b) Apply the log function to stabilize variable and analyze the variable yt=log(eps) at time t. The distribution of log(eps) is closer to normal distribution, with symmetry and kurtosis statistics

close to 0 as also shown by histogram and normal quantile plot. The Jarque-Bera test is not

significant (p-value > 0.05) indicating that hypothesis of normal distribution for log(eps) cannot

be rejected. Moments

N 48 Sum Weights 48

Mean -0.4891602 Sum Observations -23.479688

Std Deviation 0.39954911 Variance 0.15963949

Skewness 0.0058175 Kurtosis -0.398239

Uncorrected SS 18.9883843 Corrected SS 7.50305604

Coeff Variation -81.68063 Std Error Mean 0.05766995

Results of Jacque and Bera test on normality

Jarque & P-value for

skewness, kurtosis, Bera Jarque &

Obs xlog xlog statistic Bera test

1 .005817502 -0.39824 0.31746 0.85323

c) Is there any evidence of a seasonal effect in the log(eps) data? Explain your answer. The time plot for Costco quarterly earnings and for log(earnings) show quarterly seasonality, with

higher earnings in the first and third quarters and lower earnings in the second and fourth

quarters. We can observe a slight dip in earnings in 2009 but the recent economic crisis didn’t

seem to affect Costco earnings much, as earnings increased shortly after the dip. The log(eps) time

plot shows that variability is more constant over time.

eps

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

date

01/01/2000 01/01/2002 01/01/2004 01/01/2006 01/01/2008 01/01/2010 01/01/2012 01/01/2014

Costco earnings

d) Compute the Dickey-Fuller test for the quarterly earnings for p=1, 3, and 5. Write down the test hypotheses, and the test statistics and analyze the results of the test

Dickey Fuller Test is expressed as H0: φ1=1 vs Ha: φ1<>1 where the null hypothesis

indicates unit-root non-stationarity.

The Dickey-Fuller test shows that the earnings time series is unit-root non-stationary. The test p-

values for lags 1, 3 and 5 are all larger than 5, therefore the null hypothesis of non stationarity

cannot be rejected (large p-values> 0.1).

The ARIMA Procedure

Augmented Dickey-Fuller Unit Root Tests

Type Lags Rho Pr < Rho Tau Pr < Tau F Pr > F

Zero Mean 1 1.1538 0.9282 1.37 0.9554

3 1.2386 0.9381 4.43 0.9999

5 1.1075 0.9213 1.91 0.9852

Single Mean 1 -0.7726 0.9018 -0.28 0.9193 1.22 0.7642

3 1.8677 0.9962 2.06 0.9998 10.02 0.0010

5 0.4869 0.9736 0.24 0.9722 1.86 0.6076

Trend 1 -20.9931 0.0306 -2.77 0.2161 4.22 0.3620

3 0.1628 0.9958 0.06 0.9958 2.31 0.7221

5 93.4268 0.9999 -2.07 0.5472 2.53 0.6813

e) Take the first difference of log(eps) series, and analyze its autocorrelation function. Does the differenced time series show strong seasonal patterns? Explain.

The ACF of the first difference of yt=log(eps), say yt-yt-1 shows larger values at lags multiple of

four, which is consistent with the seasonal behavior seen in the time plot. The TS is not stationary

since it shows seasonality.

xlog

-1.3

-1.2

-1.1

-1.0

-0.9

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

0.1

0.2

0.3

0.4

date

01/01/2000 01/01/2002 01/01/2004 01/01/2006 01/01/2008 01/01/2010 01/01/2012 01/01/2014

Time plot of Costco log(EPS)

f) If you de-trend and de-seasonalize the Costco log(eps) data, do you obtain a stationary time series? Use appropriate statistical techniques and graphs to analyze stationarity.

The de-trended and de-seasonalized time series is stationary – as indicated by the Dickey-Fuller

test (p-values are small, and therefore hypothesis of non stationarity can be rejected), and small

ACF values. After differencing the time series, we can assume a zero mean model for the data.

Note that deseasonalized TS have zero mean (use”nc” type in R)

The ARIMA Procedure

Augmented Dickey-Fuller Unit Root Tests

Type Lags Rho Pr < Rho Tau Pr < Tau F Pr > F

Zero Mean 1 -39.9069 <.0001 -4.17 0.0001

3 85.9684 0.9999 -5.29 <.0001

5 28.6760 0.9999 -4.61 <.0001

Single Mean 1 -39.7843 0.0003 -4.10 0.0024 8.49 0.0010

3 85.8868 0.9999 -5.19 0.0002 13.58 0.0010

5 28.7344 0.9999 -4.51 0.0009 10.30 0.0010

Trend 1 -39.6292 <.0001 -4.04 0.0149 8.22 0.0184

3 86.3878 0.9999 -5.10 0.0009 13.11 0.0010

5 28.3836 0.9999 -4.48 0.0053 10.04 0.0010

g) Apply an airline model SARIMA(0,1,1)(0,1,1)4 to estimate the log(eps).

a. Write down the estimated expression of the fitted model (1-B)(1-B

4 )Xt = (1-0.339B)(1-0.669B

4 )at (sas output)

The ARIMA Procedure

Maximum Likelihood Estimation

Standard Approx

Parameter Estimate Error t Value Pr > |t| Lag

MA1,1 0.33862 0.14912 2.27 0.0232 1

MA2,1 0.66864 0.13492 4.96 <.0001 4

Variance Estimate 0.012082

Std Error Estimate 0.109917

AIC -63.4007

SBC -59.8783

Number of Residuals 43

b. Are all model coefficients significant? All coefficients of the SARIMA model are significant (p-value < 0.05). .

c. Analyze the residuals to check the model validity The autocorrelation analysis of the residuals show that the ACF at lag 5 is significant larger

than zero. This is also shown by the Ljung Box test that has p-value > 0.05 for lags up to 6.

These results indicate that there is some source of serial correlation that is still unexplained

by the model. Since it is only one large ACF value, this is not indicative of a seasonal trend.

The solution is to add a lag 5 parameter in the MA component of the SARIMA model.

Autocorrelation Check of Residuals

To Chi- Pr >

Lag Square DF ChiSq --------------------Autocorrelations--------------------

6 9.44 4 0.0509 -0.013 0.104 -0.198 -0.129 -0.346 -0.049

12 12.75 10 0.2380 0.002 0.184 0.019 0.153 -0.000 0.023

18 20.68 16 0.1911 -0.185 0.066 -0.179 -0.065 -0.044 0.189

24 21.88 22 0.4673 -0.014 -0.055 -0.066 0.054 -0.053 -0.009

h) The multiplicative airline model has some limitations, how can you modify it to be a better fit for the data? (If you are not sure how to do this, post a message on the discussion forum, and start a

discussion with other students!)

a. Write down the estimated expression of the fitted model (1-B)(1-B

4 )Xt =(1-0.7235B

4 )(1-0.4103B-0.3726B

5 )at (sas output)

The ARIMA Procedure

Maximum Likelihood Estimation

Standard Approx

Parameter Estimate Error t Value Pr > |t| Lag

MA1,1 0.41029 0.16632 2.47 0.0136 1

MA1,2 0.37262 0.16359 2.28 0.0227 5

MA2,1 0.72350 0.13731 5.27 <.0001 4

Variance Estimate 0.010536

Std Error Estimate 0.102645

AIC -66.3638

SBC -61.0802

Number of Residuals 43

b. Are all model coefficients significant? All coefficients are significantly different from zero

c. Analyze the residuals to check the model validity The residuals are white noise. Ljung Box test for lags up to 24 is significant, indicating that

the model is able to capture the time dependencies in the data.

To Chi- Pr >

Lag Square DF ChiSq --------------------Autocorrelations--------------------

6 3.77 3 0.2870 -0.059 0.091 -0.202 -0.134 -0.078 0.028

12 5.94 9 0.7457 0.026 0.140 -0.069 0.078 -0.055 0.062

18 13.34 15 0.5759 -0.109 0.071 -0.205 -0.086 -0.066 0.182

24 14.27 21 0.8577 0.005 -0.042 -0.079 0.018 -0.033 -0.032

i) Which model is a better fit for the data? Explain what evidence you use to select the best model.

The SARIMA(0,1,0)(0,1,5)4 provides a better description of the dynamic behavior of the Costco

log earnings time series. The coefficients are significantly different from zero, residuals are white

noise and normally distributed. Also the second model has the lowest SBC/BIC values.

j) Compute forecasts for Costco quarterly EPS for the next 4 quarters using the selected model in h). The table below reports forecasts in terms of EPS values. The forecasts of log(eps) computed by

SARIMA model were transformed to EPS original scale using the exp() function.

Obs date OBSERVED FORECAST STD L95 U95 RESIDUAL x

49 10/01/2013 . 0.97486 0.10264 0.79302 1.18584 . .

50 01/01/2014 . 1.15743 0.11916 0.90987 1.45159 . .

51 04/01/2014 . 1.00514 0.13366 0.76662 1.29454 . .

52 07/01/2014 . 1.51769 0.14672 1.12620 2.00170 . .

R code: # Analysis of quarterly house starts (single family homes) in the US

from 1974

# to 2012. Seasonality is expected due to lower number of new

constructions

# in the winter months.

#1. LOAD LIBRARIES

library(tseries)

library(fBasics)

library(forecast)

#2. IMPORT DATA

# Load data with no variable names into the data frame "da"

myd=read.table("costco.csv",header=T, sep=',')

# creates time series object

x=myd$eps

xts=ts(myd$eps,start=c(2001,4), frequency=4,)

# create time plot

plot(xts,type='l')

# TRANSFORMATIONS

# natural log

y=log(x)

yts=log(xts)

plot(yts,type='l')

#add circles on data points in the time plot

points(yts)

# ACF ANALYSIS

# two plots per page

par(mfcol=c(1,1))

# ACF plot for log starts

# as.vector() transforms time series to regular numerical vector

acf(as.vector(y),lag.max=24, main="ACF of log starts")

# APPLYING DIFFERENCING TO DATA

# compute regular differences

dy=diff(y)

# create acf plot

acf(as.vector(dy),lag.max=16, main="ACF of DX log starts")

# compute seasonal difference for quarterly data (s=4)

sdy=diff(dy,4)

# create acf plot

acf(as.vector(sdy),lag.max=16, main="ACF of DSDX log starts")

pacf(sdy)

#FIT SARIMA MODEL

library(forecast)

# try automated order selection

# since TS is stationary, stationary = FALSE;

auto.arima(y, trace=T,stationary = FALSE)

# fit multiplicative seasonal model ARIMA(1,0,0)×(1,1,0)[4]

m1=arima(y, order=c(1,0,0),seasonal=list(order=c(1,1,0),period=4),

method="ML")

source("arimaPrint.R")

# use arimaPrint() function in arimaPrint.R file

# save arimaPrint.R file in work directory

arimaPrint(m1)

acf(m1$residuals)

qqnorm(m1$residuals)

qqline(m1$residuals)

hist(m1$residuals)

# Since seasonal AR_1 coefficient is not significant; remove it from

model;

# fit airline seasonal model ARIMA(0,1,1)×(0,1,1)4

m1=arima(y, order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4),

method="ML")

m1

# fit airline seasonal model ARIMA(0,1,5)×(0,1,1)4 with only lags 1

and 5

m1=arima(y, order=c(0,1,5),seasonal=list(order=c(0,1,1),period=4),

method="ML", fixed=c(NA, 0,0,0,NA,NA))

m1

# use arimaPrint() function in arimaPrint.R file

arimaPrint(m1,c(1,5,6) )

# One plot per page

par(mfcol=c(1,1))

acf(m1$residuals)

# "approximated" ljung box test on residuals

Box.test(m1$residuals, 6, "Ljung-Box",fitdf=3 )

Box.test(m1$residuals, 12, "Ljung-Box", fitdf=3)

# compute predictions for log data for up to 8-step aheads

f1=forecast(m1,h=8)

names(f1)

plot(f1)

#predictions in original scale;

exp(f1$mean)

# creates forecast plot for original variable

# join data with forecasts

s1=c(x,exp(f1$mean))

#s1=ts(s1, frequency=4,start=c(1974,1))

# computes lower limit for 95% interval

lcl=c(x,exp(f1$lower[,2]))

# computes upper limit for 95% interval

ucl=c(x,exp(f1$upper[,2]))

# Forecast plot

plot(s1,type='l',col='red')

# points(s2011)

lines(1:length(s1),ucl,lty=2)

lines(1:length(s1),lcl,lty=2)

SAS code title1 'Costco earnings';

proc import datafile="costco.csv" out=myd replace;

getnames=yes;

delimiter=',';

run;

proc print;

run;

data myd;

set myd;

xlog=log(eps);

run;

*analyze distribution of log(eps);

proc univariate data=myd;

var xlog;

histogram/normal;

qqplot / normal(mu=est sigma=est);

output out=stats kurtosis=kurtosis skewness=skewness N=ntot;

run;

/* steps to compute Jarque-Bera tests*/

/* and print out tests results to the OUTPUT window */

data computation;

set stats;

label pv_kur = "P-value for kurtosis test";

*skewness statistic;

skew_test = skewness/sqrt(6/Ntot);

*kurtosis statistic;

kurt_test = kurtosis/sqrt(24/Ntot);

* jarque-bera statistic;

jb = skew_test*skew_test+kurt_test*kurt_test;

* JB test p-values;

pv_jb = 1-cdf('CHISQUARE', jb,2);

label pv_jb = "P-value for Jarque & Bera test"

jb = "Jarque & Bera statistic";

/* Print out results of tests*/

Title " Results of Jacque and Bera test on normality";

proc print data= computation label;

var skewness kurtosis jb pv_jb;

run;

*plot timeplots for all transformed variables;

symbol1 i=join v=dot;

proc gplot data=myd;

title1 "Time plot of Costco EPS";

plot eps * date = 1 ;

title1 "Time plot of Costco log(EPS)";

plot xlog * date = 1 ;

run;

proc arima data=myd;

identify var=eps nlag=24 stationarity=(adf=(1 3 5));

identify var=xlog(1) nlag=24;

identify var=xlog(1,4) nlag=20 stationarity=(adf=(1 3 5))

minic;

identify var=xlog(1,4) nlag=20 esacf;

run;

*airline model SARIMA(0,1,1)(0,1,1)[4];

estimate q=(1)(4) noconstant method=ml;

* modified airline model SARIMA(0,1,1)(0,1,5)[4];

estimate q=(1 5)(4) noconstant method=ml;

outlier maxnum=3 alpha=0.01;

run;

forecast out=b lead=12 id=date interval=quarter;

quit;

/*

The following statements retransform the forecast values to get

forecasts in the original scales.

*/

data c;

set b;

x = exp( xlog );

forecast = exp( forecast + std*std/2 );

l95 = exp( l95 );

u95 = exp( u95 );

run;

proc print data=c;

run;

Problem 3 [12 points] (review Week 5 notes)

Global temperatures are rising. Analyze the time series of Land-Ocean Temperature Index in degrees

Celsius published at the NASA-GISS (Goddard Institute for Space Studies) and saved in dataset

globtemp.csv. The temperature index measures temperature anomalies that are computed relative

to the base period 1951-1980. The reason to work with anomalies, rather than absolute

temperature is that absolute temperature varies markedly in short distances, while monthly or

annual temperature anomalies are representative of a much larger region (from GISS website). a) Plot the temperature index time series and its ACFs (20 lags). Analyze trends and patterns

shown by the data.

Time plot shows that global temperatures has steadily increased since 1950.The ACF

plot shows that Global Temperature Index TS is non stationary (ACF function decays

very slowly to zero, which is not surprising given the upward time trend shown by the

time plot)

b) Analyze if the series is stationary using both the ACF function and the Dickey Fuller test The DF tests computed for AR(3) and AR(5) models show that the null hypothesis of unit-root

stationarity cannot be rejected (p-values are > 0.1). (NOTE: for this dataset, it would be

acceptable to use the DF test for either model with intercept (constant) but no time trend, or

for a model with intercept (constant) and a time trend. Results are equivalent)

> # Dickey Fuller test for time series with time trend

> library(fUnitRoots)

> adfTest(index, lags=3, type=c("ct"))

Title:

Augmented Dickey-Fuller Test

Test Results:

PARAMETER:

Lag Order: 3

STATISTIC:

Dickey-Fuller: -2.4227

P VALUE:

0.4007

> adfTest(index, lags=5, type=c("ct"))

Title:

Augmented Dickey-Fuller Test

Test Results:

PARAMETER:

Lag Order: 5

STATISTIC:

Dickey-Fuller: -1.9472

P VALUE:

0.5984

In SAS

Augmented Dickey-Fuller Unit Root Tests

Type Lags Rho Pr < Rho Tau Pr < Tau F Pr > F

Zero Mean 1 -3.1286 0.2229 -0.97 0.2967

3 0.3095 0.7561 0.16 0.7308

5 0.6287 0.8363 0.36 0.7864

Single Mean 1 -2.9490 0.6593 -0.91 0.7827 0.74 0.8811

3 0.5894 0.9790 0.31 0.9780 1.13 0.7831

5 1.0265 0.9887 0.63 0.9900 1.84 0.6012

Trend 1 -26.5191 0.0132 -3.75 0.0225 7.46 0.0207

3 -11.2587 0.3399 -2.42 0.3659 4.15 0.3489

5 -8.1422 0.5607 -1.95 0.6245 3.26 0.5268

c) Use order selection methods such as EACF or BIC to identify the order of the “best” AR/MA/ARMA model for the global temperature index time series. Note that since the TS is

not stationary, you should use an ARIMA(p,1,d) model.

Since TS is unit-root non stationary, only ARIMA(p,1,q) models will be evaluated.

In R: The EACF method suggests an ARIMA(1,1,2). The auto.arima() function selects the

ARIMA(0,1,1) model that exhibits the lowest BIC value.

In SAS: the EACF method suggests an ARIMA(0,1,1) model and the BIC model selects the

ARIMA(3,1,0)

> #EACF method > library(TSA) > #use numeric object not ts object. coredata() function retrieves numeric vector > d1_ts=diff(indexts) > eacf(coredata(d1_ts)) AR/MA 0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 x o o x o o o o x o o o x o 1 x o o o o o o o o o o o o o

2 x o x o o o o o o o o o x o 3 o x x o o o o o o o o o o o 4 o x x x o o o o o o o o o o 5 x o x x x o o o o o o o o o 6 x o o o x o o o o o o o o o 7 x o o o x x o o o o o o o o

> # BIC criterion > # use numeric object not ts object > auto.arima(index, max.p = 5, max.q = 5, stationary = F, ic =c("bic"), trace=TRUE, allowdrift=F) ARIMA(2,1,2) : 982.1547 ARIMA(0,1,0) : 990.837 ARIMA(1,1,0) : 985.1276 ARIMA(0,1,1) : 972.8737 ARIMA(1,1,1) : 976.1241 ARIMA(0,1,2) : 973.7256 ARIMA(1,1,2) : 979.3824 Best model: ARIMA(0,1,1) Series: index ARIMA(0,1,1) Coefficients: ma1 -0.5398 s.e. 0.0838 sigma^2 estimated as 89.26: log likelihood=-483.91 AIC=971.82 AICc=971.92 BIC=977.59

In SAS ESACF Probability Values

Lags MA 0 MA 1 MA 2 MA 3 MA 4 MA 5

AR 0 0.0011 0.0635 0.4021 0.0552 0.1309 0.8882

AR 1 <.0001 0.7810 0.3314 0.3416 0.1542 0.5998

AR 2 <.0001 0.6174 0.0018 0.7358 0.1896 0.8318

AR 3 0.4231 0.0002 0.0001 0.1756 0.6531 0.7922

AR 4 0.2412 <.0001 0.0018 0.0597 0.3335 0.6631

AR 5 <.0001 0.4323 0.0527 0.0004 0.0120 0.9065

The ARIMA Procedure

ARMA(p+d,q)

Tentative

Order

Selection

Tests

---ESACF---

p+d q

0 1

3 3

(5% Significance Level)

Name of Variable = temp

Minimum Information Criterion

Lags MA 0 MA 1 MA 2 MA 3 MA 4 MA 5

AR 0 4.588787 4.487776 4.468956 4.495152 4.505573 4.533958

AR 1 4.541124 4.523001 4.480694 4.48462 4.503237 4.540091

AR 2 4.499424 4.522147 4.516266 4.503816 4.537365 4.574281

AR 3 4.441992 4.477332 4.494755 4.527052 4.558184 4.586585

AR 4 4.477766 4.514051 4.522759 4.556779 4.577564 4.612352

AR 5 4.484678 4.521373 4.55314 4.590039 4.610067 4.637321

Error series model: AR(10)

Minimum Table Value: BIC(3,0) = 4.441992

d) Fit the selected ARIMA model, and write down the result. The ARIMA(1,1,2) model selected by the EACF method in R has non-significant parameters.

The ARIMA(3,1,0) selected by BIC method in SAS is a possible model but has higher BIC

value (details of analysis are omitted). The best model is the ARIMA(0,1,1) selected by the BIC

criterion in R and EACF in SAS.

The fitted ARIMA(0,1,1) model for the global temperature index Xt can be written as:

(1-B)Xt = (1-0.5398B)at

> m2= Arima(index, order=c(0,1,1), method='ML') > m2 Series: index ARIMA(0,1,1) Coefficients: ma1 -0.5398 s.e. 0.0838 sigma^2 estimated as 89.26: log likelihood=-483.91 AIC=971.82 AICc=971.92 BIC=977.59 > arimaPrint(m2) [1] "T-test table for model coefficients" estimates std error ma1 "-0.5397611" "0.08376795" p-values ma1 "1.167279e-10"

e) Use residual analysis to check the model. Ljung Box test on residuals shows that hypothesis of

independence cannot be rejected. So we can assume that

residuals are white noise.

The normal quantile plot shows that distribution of residuals

is close to the normal distribution.

> Box.test(m2$residuals,lag=6,type='Ljung-Box', fitdf=1)

Box-Ljung test

data: m2$residuals

X-squared = 9.5895, df = 5,

p-value = 0.08774

> Box.test(m2$residuals,lag=12,type='Ljung-Box', fitdf=1)

Box-Ljung test

data: m2$residuals

X-squared = 16.3747, df = 11,

p-value = 0.1278

f) Compute 1-step to 4-step ahead forecasts of the fitted model with origin at the end of the data, i.e. 2012. Write down the forecasts and their standard errors.

> f=forecast.Arima(m2, h=4) > plot(f, include=50) > f Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 134 58.1793 46.07184 70.28676 39.66253 76.69606 135 58.1793 44.85108 71.50751 37.79555 78.56305 136 58.1793 43.73311 72.62548 36.08577 80.27283 137 58.1793 42.69566 73.66293 34.49912 81.85947

g) EXTRA CREDIT [2pts]: apply backtesting procedure to compare the following two models: ARIMA(0,1,1) and ARIMA(0,1,2). Use 113

cases (85% of data) for training and 20 cases

(or 15%) for testing.

The procedure shows that the RMSE value of

the ARIMA(0,1,1) is 9.680. While the RMSE

value for the ARIMA(0,1,2) is 9.688. The two

values are very close showing that the two

models are equivalent in terms of predictive

power. However the ARIMA(0,1,1) is slightly

better.

> backtest(m2,index,113,1) [1] "RMSE of out-of-sample forecasts" [1] 9.680175 [1] "Mean absolute error of out-of-sample forecasts" [1] 7.925004 > m3=arima(index, order=c(0,1,2)) > backtest(m3,index,113,1) [1] "RMSE of out-of-sample forecasts" [1] 9.688697 [1] "Mean absolute error of out-of-sample forecasts" [1] 7.932232 > m3 Series: x ARIMA(0,1,2) Coefficients: ma1 ma2 -0.4287 -0.1591 s.e. 0.0843 0.0764

sigma^2 estimated as 86.63: log likelihood=-481.94 AIC=967.88 AICc=968.07 BIC=976.53

In SAS: backtest needs to be modified to fit ARIMA(p,1,q) models by changing the IDENTIFY

statement in the PROC ARIMA of the BACKTEST macro as follows identify var=&VAR(1) noprint;

/*--------- backtesting --------------*/

%backtest(trainpct=85, dataset=myd, var=temp, q=1, interval=year);

%backtest(trainpct=85, dataset=myd, var=temp, q=2, interval=year);

Backtest results for myd

Model: q=(1)

Obs _TYPE_ _FREQ_ mafe msfe rmsfe

1 0 20 7.71428 87.0438 9.32973

Backtest results for myd

Model: q=2

Obs _TYPE_ _FREQ_ mafe msfe rmsfe

1 0 20 7.60769 85.3367 9.23779

R code

# Analysis of global temperature index

#1. LOAD LIBRARIES

library(tseries)

library(fBasics)

library(forecast)

#Create a time series

library(zoo)

data=read.table("globtemp.csv",header=T, sep=',')

index=data[,2]

indexts=zoo(data[,2],as.Date(as.character(data[,1]),format="%Y"))

#plot ts and acf test

plot(indexts, type='l', xlab='time', ylab='temperature index ')

#Analysis of distribution

par(mfcol=c(1,1))

hist(index, xlab="Distribution of GT Index", prob=TRUE,

main="Histogram")

# add approximating normal density curve

xfit<-seq(min(index),max(index),length=40)

yfit<-dnorm(xfit,mean=mean(index),sd=sd(index))

lines(xfit, yfit, col="blue", lwd=2)

# CREATE NORMAL PROBABILITY PLOT

qqnorm(index)

qqline(index, col = 2)

acf(coredata(indexts), plot=T, lag=20)

# Dickey Fuller test

library(fUnitRoots)

adfTest(index, lags=3, type=c("ct"))

adfTest(index, lags=5, type=c("ct"))

adfTest(index, lags=3, type=c("c"))

adfTest(index, lags=5, type=c("c"))

#Ljung - Box test

Box.test(index,lag=6,type='Ljung')

Box.test(index,lag=12,type='Ljung')

#EACF method

library(TSA)

#use numeric object not ts object. coredata() function retrieves

numeric vector

d1_ts=diff(indexts)

eacf(coredata(d1_ts))

# BIC criterion

# use numeric object not ts object

auto.arima(index, max.p = 5, max.q = 5, stationary = F, ic

=c("bic"), trace=TRUE, allowdrift=F)

# 13. MODEL FITTING

# MODEL 1:

#EACF suggests ARIMA(1,1,2) p=1 and m=2

m1= arima(index, order=c(1,1,2), method='ML')

m1

source("arimaPrint.R")

arimaPrint(m1)

# RESIDUAL ANALYSIS

Box.test(m1$residuals,lag=6,type='Ljung-Box', fitdf=1)

Box.test(m1$residuals,lag=12,type='Ljung-Box', fitdf=1)

acf(m1$residuals)

qqnorm(m1$residuals)

qqline(m1$residuals)

hist(m1$residuals)

# MODEL 2 SELECTED BY AUTOARIMA

# FIT AN ARIMA(0,1,1) MODEL

#

m2= Arima(index, order=c(0,1,1), method='ML')

m2

arimaPrint(m2)

# RESIDUAL ANALYSIS

Box.test(m2$residuals,lag=6,type='Ljung-Box', fitdf=1)

Box.test(m2$residuals,lag=12,type='Ljung-Box', fitdf=1)

acf(m2$residuals)

#Analysis of distribution

par(mfcol=c(1,1))

hist(m2$residuals, xlab="Distribution of residuals", prob=TRUE,

main="Histogram")

# add approximating normal density curve

xfit<-seq(min(m2$residuals),max(m2$residuals),length=40)

yfit<-dnorm(xfit,mean=mean(m2$residuals),sd=sd(m2$residuals))

lines(xfit, yfit, col="blue", lwd=2)

# CREATE NORMAL PROBABILITY PLOT

qqnorm(m2$residuals)

qqline(m2$residuals, col = 2)

#forecast and plot

pred=predict(m2,n.ahead=4, se.fit=T)

pred

forecast.Arima(m2, h=4)

f=forecast.Arima(m2, h=4)

plot(f, include=50)

#Backtesting

source("backtest.R")

backtest(m2,index,113,1)

m3=arima(index, order=c(0,1,2))

backtest(m3,index,113,1)

SAS code: %macro backtest(TRAINPCT=80, DATASET=, VAR=, P=, Q=, D=, DATE=date,

INTERVAL=month);

* -------------------------------------------------------------------------

;

* T I M E S E R I E S B A C K T E S T M A C R O

* -------------------------------------------------------------------------

;

* DePaul CSC425, Winter 2013, Dr. Raffaella Settimi

* Macro written by Bill Qualls, First Analytics

* -------------------------------------------------------------------------

;

* E X P L A N A T I O N O F P A R A M E T E R S

* (Order of variables is insignificant)

* -------------------------------------------------------------------------

;

* TRAINPCT .. Percent of dataset to be used for training.

* So, (100 - TRAINPCT)% will be used for evaluation.

* Example: TRAINPCT=80

* DATASET .. Time series dataset. Libname optional, defaults to Work.

* Example: DATASET=Work.Unemp

* VAR .. Name of time series variable.

* Example: VAR=ratechg

* P .. Specified for AR models. Omit otherwise.

* Example: P=(1 3 6)

* Q .. Specified for MA models. Omit otherwise.

* Example: Q=(1 3 6)

* DATE .. Name of date variable. Defaults to date.

* Example: DATE=date

* INTERVAL .. Date interval. Defaults to month.

* Example: INTERVAL=day

* -------------------------------------------------------------------------

;

* S A M P L E U S A G E

* %backtest(trainpct=80, dataset=work.unemp, var=ratedif, p=(1),

interval=day);

* -------------------------------------------------------------------------

;

* How many records are in the dataset? ;

data _null_;

call symput('NRECS', trim(left(nrecs)));

set &DATASET nobs=nrecs;

stop;

run;

* Determine which ones are exclusively for training based on TRAINPCT ;

%let SIZE_OF_ROLLING_WINDOW = %sysfunc(round(&NRECS * &TRAINPCT / 100));

* create a working copy of dataset with observation number ;

* as a variable for use in a where clause with proc arima. ;

* also add a placeholder for the predicted value. ;

data Work._MY_COPY_ (keep = &VAR &DATE _OBS_ _PRED_);

set &DATASET;

_OBS_ = _N_;

_PRED_ = .;

run;

* turn off log -- too lengthy ;

filename junk dummy;

proc printto log=junk print=junk;

run;

* Will build the model once for each record used in evaluation. ;

* Each time I will predict one record forward. ;

%let MODELS_TO_BE_BUILT = %sysevalf(&NRECS - &SIZE_OF_ROLLING_WINDOW);

%do i = 1 %to &MODELS_TO_BE_BUILT;

* Model using SIZE_OF_ROLLING_WINDOW records, and make one prediction ;

proc arima data=Work._MY_COPY_ plots=none;

where _OBS_ ge &i

and _OBS_ le (&i + &SIZE_OF_ROLLING_WINDOW - 1);

identify var=&VAR(1) noprint;

estimate

%if (&P ne ) %then %do; p=&P %end;

%if (&Q ne ) %then %do; q=&Q %end;

method=ml noprint;

forecast lead=1 id=&DATE interval=&INTERVAL out=Work._MY_RESULTS_

noprint;

run;

* get the predicted value (in the last record) as a macro variable ;

data _null_;

p = nrecs;

set Work._MY_RESULTS_ point=p nobs=nrecs;

call symput("FORECAST", forecast);

stop;

run;

* move that prediction to its place in the output file ;

proc sql noprint;

update Work._MY_COPY_

set _PRED_ = &FORECAST

where _OBS_ = &i + &SIZE_OF_ROLLING_WINDOW;

quit;

run;

* show progress so far ;

%if (%sysfunc(mod(&i, 5)) = 0) %then %do;

* print on;

proc printto log=log print=print;

run;

%put Finished &i iterations;

* print off again;

proc printto log=junk print=junk;

run;

%end;

%end;

* turn print back on ;

proc printto log=log print=print;

run;

* calculate prediction error;

data Work._MY_COPY_;

set Work._MY_COPY_;

Predicted_Error = (&VAR - _PRED_) ;

Predicted_Error_Squared = (&VAR - _PRED_)**2;

absresidual = abs(Predicted_Error);

run;

* compute and report the mean square forecast error;

%if (&P eq ) %then %let PP = ; %else %let PP = p=&P;

%if (&Q eq ) %then %let QQ = ; %else %let QQ = q=&Q;

title "Backtest results for &DATASET";

title2 "Model: &PP &QQ";

proc summary data=Work._MY_COPY_;

where _OBS_ > &SIZE_OF_ROLLING_WINDOW;

var Predicted_Error absresidual;

output out=outm mean(absresidual)=mafe mean(Predicted_Error_Squared)=msfe;

run;

data outm;

set outm;

rmsfe=sqrt(msfe);

run;

proc print data=outm;

run;

/*

proc means data=Work._MY_COPY_ N MEAN STDDEV CLM T PROBT;

where _OBS_ > &SIZE_OF_ROLLING_WINDOW;

var Predicted_Error_Squared;

run;

title;

*/

%mend backtest;

/*-------------------------------------------------*/

proc import datafile="globtemp.csv" out=myd replace;

delimiter=',';

getnames=yes;

run;

proc print data=myd;

run;

* produced time plots;

goptions reset=global;

symbol interpol = join value=star;

proc gplot data=myd;

plot temp*year / haxis = 1880 to 2012 by 10;

run;

quit;

goptions reset = global;

proc univariate data=myd;

var temp;

histogram / normal;

qqplot / normal(mu=est sigma=est);

run;

* compute Dickey-Fuller unit-root tests;

proc arima data=myd;

identify var = temp stationarity = (adf=(1 3 5));

identify var = temp(1) stationarity = (adf=(1 3 5));

run;

proc arima data=myd;

identify var = temp esacf q=(0:5) p=(0:5) nlag =20;

identify var = temp minic q=(0:5) p=(0:5) nlag =20;

run; quit;

/* fit models from step above */

* fit ARIMA(p,1,q) models for unemployment rates;

* note that model intercept is set =0 using the noconstant option;

proc arima data=myd;

identify var = temp(1) noprint;

estimate p=0 q=1 method=ml plot noconstant;

estimate p=3 q=0 method=ml plot noconstant;

run;

quit;

proc arima data=myd;

identify var=temp(1) noprint;

estimate p=0 q=1 method=ml plot noconstant;

forecast back=0 lead=4 printall out=result id=year interval=day;

run;

quit;

/* Analysis of residuals*/

/* Residual plot: residual vs time*/

goptions reset=global;

symbol i=join;

proc gplot data=result;

plot residual*year / haxis = 1880 to 2012 by 10;

run;

quit;

/*analysis of distribution of residuals, to check

for normality */

goptions reset=global;

proc univariate data=result;

var residual;

histogram/normal (mu=est sigma=est);

qqplot/normal(mu=est sigma=est);

run;

goptions reset=global;

symbol1 interpol =join value=star color=red w=2;

symbol2 interpol = join color=black;

symbol3 interpol = join value=none color = blue;

symbol4 interpol = join value=none color = blue;

legend1 label=none

shape=symbol(4,2)

position=(top center inside)

mode=share;

proc gplot data=result;

plot forecast*year temp*year u95*year l95*year / overlay

haxis = 1880 to 2022 by 10 legend=legend1;

run;

/*--------- backtesting --------------*/

%backtest(trainpct=85, dataset=myd, var=temp, q=(1), interval=year);

%backtest(trainpct=85, dataset=myd, var=temp, q= 2, interval=year);