setwd("~/teaching/time-series/ts2015spring/data")

###----------------------------------------------------------------------------
### ARIMA model
###----------------------------------------------------------------------------

da=read.table("m-ibm3dx2608.txt",header=T)
sibm=da[,2] # Get the IBM simple returns
Box.test(sibm,lag=5,type='Ljung') # Ljung-Box statistic Q(5) Box-Ljung test

## Example 2.1
gnp=scan(file='dgnp82.txt')
gnp1=ts(gnp,frequency=4,start=c(1947,2))
plot(gnp1)
points(gnp1,pch='*')
m1=ar(gnp,method="mle") # Find the AR order
m1$order # An AR(3) is selected based on AIC
m2=arima(gnp,order=c(3,0,0)) # Estimation
m2

## In R, ‘‘intercept'' denotes the mean of the series.
## Therefore, the constant term is obtained below:
(1-.348-.1793+.1423)*0.0077

p1=c(1,-m2$coef[1:3]) # Characteristic equation
roots=polyroot(p1) # Find solutions


###----------------------------------------------------------------------------
### Dickey-Fuller Test (Example 2.2)
###----------------------------------------------------------------------------
library(fUnitRoots)
da=read.table("q-gdp4708.txt",header=T)
plot(da[, c(1, 4)],  type = "l")

gdp=log(da[,4])
plot(da[,1], log(da[, 4]),  type = "l")

plot(diff(da[, 4], 2),  type = "l")

m1=ar(diff(gdp),method='mle')

adfTest(diff(gdp),lags=10,type=c("c"))


###----------------------------------------------------------------------------
### Testing for autocorrelation
###----------------------------------------------------------------------------
library("FinTS")
da=read.table("m-intc7308.txt",header=T)
intc=log(da[,2]+1)
AutocorTest(intc,lag=12)

###----------------------------------------------------------------------------
### Seasional time series (Example 2.4)
###----------------------------------------------------------------------------
da=read.table("m-deciles08.txt",header=T)
d1=da[,2]
jan=rep(c(1,rep(0,11)),39) # Create January dummy.
plot(da[, 2], type = "l", col = "red")
plot(d1[-1]-lag(d1, 12), type = "l")
m1=lm(d1~jan)
summary(m1)

m2=arima(d1,order=c(1,0,0),seasonal=list(order=c(1,0,1),period=12))
m2
arima(d1, order = c(1, 0, 0))
tsdiag(m2,gof=36)

m2=arima(d1,order=c(1,0,0),seasonal=list(order=c(1,0,1), period=12),include.mean=F)
m2