#getting the data for WI, US, IL, MI, MN, IO #Using Quandl data, great little site wi<-read.csv('http://www.quandl.com/api/v1/datasets/FRED/WIUR.csv?&auth_token=gigXwpxd6Ex91cjgz1B7&trim_start=1976-01-01&trim_end=2013-04-01&sort_order=desc', colClasses=c('Date'='Date')) us<-read.csv('http://www.quandl.com/api/v1/datasets/FRED/UNRATE.csv?&auth_token=gigXwpxd6Ex91cjgz1B7&trim_start=1976-01-01&trim_end=2013-04-01&sort_order=desc', colClasses=c('Date'='Date')) il<-read.csv('http://www.quandl.com/api/v1/datasets/FRED/ILUR.csv?&auth_token=gigXwpxd6Ex91cjgz1B7&trim_start=1976-01-01&trim_end=2013-04-01&sort_order=desc', colClasses=c('Date'='Date')) mi<-read.csv('http://www.quandl.com/api/v1/datasets/FRED/MIUR.csv?&auth_token=gigXwpxd6Ex91cjgz1B7&trim_start=1976-01-01&trim_end=2013-04-01&sort_order=desc', colClasses=c('Date'='Date')) mn<-read.csv('http://www.quandl.com/api/v1/datasets/FRED/MNUR.csv?&auth_token=gigXwpxd6Ex91cjgz1B7&trim_start=1976-01-01&trim_end=2013-04-01&sort_order=desc', colClasses=c('Date'='Date')) io<-read.csv('http://www.quandl.com/api/v1/datasets/FRED/IAUR.csv?&auth_token=gigXwpxd6Ex91cjgz1B7&trim_start=1976-01-01&trim_end=2013-04-01&sort_order=desc', colClasses=c('Date'='Date')) #merging the data into one dataframe #I started with WI so that I could isolate the rest of the variables better unemp<-merge(wi, io, by='Date'); colnames(unemp)<-c('Date', 'wi', 'io') unemp<-merge(unemp, mi, by='Date'); colnames(unemp)<-c('Date', 'wi', 'io', 'mi') unemp<-merge(unemp, mn, by='Date'); colnames(unemp)<-c('Date', 'wi', 'io', 'mi', 'mn') unemp<-merge(unemp, us, by='Date'); colnames(unemp)<-c('Date', 'wi', 'io', 'mi', 'mn', 'us') unemp<-merge(unemp, il, by='Date'); colnames(unemp)<-c('Date', 'wi', 'io', 'mi', 'mn', 'us', 'il') save(unemp, file='Unemployment WI.RData') library(car) library(lmtest) library(forecast) load('Unemployment WI.RData') unemp.b<-unemp[1:436,] unemp.p<-unemp[437:448,] plot(unemp$Date, unemp$us, col='black',

main='Unemployment WI Region', type='l', xlab='Date',

ylab='Unemp. Rate', ylim=c(2, 17)) lines(unemp$Date, unemp$il, col='blue') lines(unemp$Date, unemp$io, col='dark green') lines(unemp$Date, unemp$mi, col='purple') lines(unemp$Date, unemp$mn, col='gray') lines(unemp$Date, unemp$wi, col='red') leg.txt<-c('US', 'IL', 'IO', 'MI', 'MN', "WI") leg.col<-c('black','blue', 'dark green', 'purple', 'gray', 'red') legend('topleft', leg.txt, col=leg.col, pch=15) summary(unemp) #correlation matrix cm<-cor(unemp[,-1]) cm #From the graph and correlation matrix, everything

`#is highly correlated, the lowest r value is 0.7879 `

#between Iowa and US. #plots of each of the data, then each individual,

`#everything looks really good, nice linear trends, `

*#although Michigan has quite bit of volatility *

#compared to the other variables. plot(unemp.b, pch='.') model<-lm(wi~., data=unemp.b[,-1]) step.m<-step(model, direction='backward') summary(step.m) #testing for outliers outlierTest(step.m) # After conducting the outlier test, looks like it is michigan #Sept 1980 where unemployment was 13.2 percent, point 57 is #close at 12.9 those are the two biggest, but the outlierTest #isolated row 56 dwtest(step.m) dwtest(step.m, alternative='two.sided') #yep lots of autocorrelation, so a multiple ARAIMA might be in order pred<-unemp.p[3:7] mult.f<-predict(step.m, newdata=pred, interval='prediction') a1<-accuracy(mult.f, unemp.p$wi) #because the MASE is missing I needed to add the following to compare later a1<-c(a1,'0') names(a1)<-c("ME", "RMSE", "MAE", "MPE", "MAPE", "MASE") # ARIMA auto.r<-auto.arima(unemp.b$wi, xreg=unemp.b[,3:7]) auto.r auto.f<-forecast(auto.r, xreg=unemp.p[,3:7], h=12) a2<-accuracy(auto.f, unemp.p$wi) a2 #Just remember that the 0 is just a place holder r.table<-rbind(a1, a2) row.names(r.table)<-c('Multiple Regression', 'Reg. with ARIMA Errors') r.table

At this point I combined the two accuracy tables to see which one was the best.

b.table<-rbind(a.table[,1:6], r.table) b.table<-b.table[order(b.table$MAPE),] b.table

The winner is: Regression with ARIMA Errors, for most of the analysis the simple mean like predictions dominated, even the multiple regression is considerably low on the list.

Good use of Quandl!

ReplyDeleteGetting Dirtier. Yes ARIMA is found to be helpful, but there are other things going on too. Due to a max of 4,096 I had to edit out about 20 outliers from the model

A change in variance was detected using TSAY's variance change detection test, automatically beginning at 9/1990 and WLS was used. We don't show the H-Weights here.

The model includes SOME of the causals(no leads or lags were detected) and outliers and three changes in seasonality(August and two in January). It also has an AR1 of .704. The forecasts were as follows with a 1.22% MAPE.

7

7

7

7

6.9

6.8

6.7

6.7

7

7.1

7.1

7.1

[(1-B**1)]Y(T) = +[X1(T)][(1-B**1)][(+ .0330)] io

+[X2(T)][(1-B**1)][(+ .274)] mn

+[X3(T)][(1-B**1)][(+ .0429)] us

+[X4(T)][(1-B**1)][(+ .206)] il

+[X5(T)][(1-B**1)][(+ 1.5171)] :PULSE 56 1980/ 8

+[X6(T)][(1-B**1)][(- .751)] :PULSE 144 1987/ 12

+[X7(T)][(1-B**1)][(+ .0258)] :PULSE 51 1980/ 3

+[X8(T)][(1-B**1)][(+ .0347)] :PULSE 52 1980/ 4

+[X9(T)][(1-B**1)][(- .0432)] :PULSE 205 1993/ 1

+[X12(T)[(1-B**1)][(- .114)] :PULSE 87 1983/ 3

+[X13(T)[(1-B**1)][(- .113)] :PULSE 88 1983/ 4

+[X37(T)[(1-B**1)][(- .0702)] :PULSE 190 1991/ 10

+[X38(T)[(1-B**1)][(- .0440)] :PULSE 382 2007/ 10

+[X39(T)[(1-B**1)][(+ .0536)] :PULSE 111 1985/ 3

+[X40(T)[(1-B**1)][(- .110)] :PULSE 42 1979/ 6

+[X41(T)[(1-B**1)][(+ .0757)] :SEASONAL PULSE 409 2010/ 1

+[X42(T)[(1-B**1)][(- .0212)] :SEASONAL PULSE 121 1986/ 1

+[X43(T)[(1-B**1)][(+ .0011)] :PULSE 192 1991/ 12

+ [(1- .704B** 1)]**-1 [A(T)]

Here is the Variance Change Detection :

DIAGNOSTIC CHECK #5: THE TSAY VARIANCE CONSTANCY TEST

The Critical value used for this test : .01

The minimum group or interval size was: 87

DIRECTION TIME DATE F VALUE P VALUE

(T)

DECREASING 177 1990/ 9 .335430 .0000

Since the automatic model fixup option for the variance stability test is

enabled, the program will now estimate the parameters of the model with a set

of weights that adjusts the residuals to account for the variance change(s).

Tom Reilly

www.autobox.com

I get an error on the line below:

ReplyDeletemult.f<-predict(step.m, newdata=pred, interval='prediction')

a1<-accuracy(mult.f, unemp.p$wi)

Error: $ operator is invalid for atomic vectors

How can I fix this? Thanks