### SOME TIME SERIES COMMANDS in R                      Last updated: 1/25/01
### ---------------------------------------------------------------------
### Note:  in R, all lines starting with # are comments, 
### i.e. they're ignored.


### LOAD THE TS LIBRARY:
library(ts)

#library()  ## lists all libraries; 

#library(help=ts) ## shows all functions in the ts library






### Input data from Chatfield page 248
x <- 
c(2.22, 2.23, 2.22, 2.20, 2.09, 1.97, 2.03, 1.98, 1.94, 1.79, 1.74, 1.86, 
1.78, 1.72, 1.79, 1.82, 1.89, 1.99, 1.89, 1.83, 1.71, 1.70, 1.97, 2.21,
2.36, 2.41, 2.92, 3.15, 3.26, 3.51, 3.48, 3.16, 3.01, 2.97, 2.88, 2.91,
3.45, 3.29, 3.17, 3.09, 3.02, 2.99, 2.97, 2.94, 2.84, 2.85, 2.86, 2.89,
2.93, 2.93, 2.87, 2.82, 2.63, 2.33, 2.22, 2.15, 2.28, 2.28, 2.06, 2.54,
2.29, 2.66, 3.03, 3.17, 3.83, 3.99, 4.11, 4.51, 4.66, 4.37, 4.45, 4.58,
4.58, 4.76, 4.89, 4.65, 4.51, 4.65, 4.52, 4.52, 4.57, 4.65, 4.74, 5.10,
5.00, 4.74, 4.79, 4.83, 4.80, 4.83, 4.77, 4.80, 5.38, 6.18, 6.02, 5.91,
5.66, 5.42, 5.06, 4.70, 4.73, 4.64, 4.62, 4.48, 4.43, 4.33, 4.32, 4.30,
4.26, 4.02, 4.06, 4.08, 4.09, 4.14, 4.15, 4.20, 4.30, 4.26, 4.15, 4.27,
4.69, 4.72, 4.92, 5.10, 5.20, 5.56, 6.08, 6.13, 6.09, 5.99, 5.58, 5.59,
5.42, 5.30, 5.44, 5.32, 5.21, 5.47, 5.96, 6.50, 6.48, 6.00, 5.83, 5.91,
5.98, 5.91, 5.64, 5.49, 5.43, 5.33, 5.22, 5.03, 4.74, 4.55, 4.68, 4.53,
4.67, 4.81, 4.98, 5.00, 4.94, 4.84, 4.76, 4.67, 4.51, 4.42, 4.53, 4.70,
4.75, 4.90, 5.06, 4.99, 4.96, 5.03, 5.22, 5.47, 5.45, 5.48, 5.57, 6.33,
6.67, 6.52, 6.60, 6.78, 6.79, 6.83, 6.91, 6.93, 6.65, 6.53, 6.50, 6.69,
6.58, 6.42, 6.79, 6.82, 6.76, 6.88, 7.22, 7.41, 7.27, 7.03, 7.09, 7.18,
6.69, 6.50, 6.46, 6.35, 6.31, 6.41, 6.60, 6.57, 6.59, 6.80, 7.16, 7.51,
7.52, 7.40, 7.48, 7.42, 7.53, 7.75, 7.80, 7.63, 7.51, 7.49, 7.64, 7.92,
8.10, 8.18, 8.52, 8.56, 9.00, 9.34, 9.04, 9.08, 9.14, 8.99, 8.96, 8.86,
8.79, 8.62, 8.29, 8.05, 8.00, 7.89, 7.48, 7.31, 7.42, 7.51, 7.71, 7.99)

### time-plot
##postscript("1to4.ps")
par(mfrow=c(2,2))
plot(c(1:252)/12,x,xlab="", ylab="",type="l")
mtext(s=1,l=2,"year")
mtext(s=2,l=2,"short-term govt security yield (%)")
mtext(s=3,l=1.5, "Figure 1:  Time plot of monthly security yield data")
abline(v=17.1)

### portion of the data
x1 <- x[1:(17*12)]
l1 <- length(x1)

### correlogram
y <- acf(x1, type="correlation", lag.max = 25, plot=F)
y1 <- y$acf
y2 <- c(y1)
plot(c(1,25),c(0,1),type="n",xlab="",ylab="")
lines(c(1:25),y2[2:26])
mtext(s=1,l=2,"lag (months)")
mtext(s=2,l=2,"autocorrelation")
mtext(s=3,l=1.5, "Figure 2:  Correlogram of (early portion of) raw data")


### take first differences
x2 <- diff(x1)
l2 <- length(x2)

### time-plot again
plot(c(0,l2/12),c(-1,1),type="n",xlab="",ylab="")
lines(c(1:l2)/12,x2)
mtext(s=1,l=2,"year")
mtext(s=2,l=2,"1st differences (%)")
mtext(s=3,l=1.5, "Figure 3:  Time plot of first differences")


### correlogram again
y <- acf(x2, type="correlation", lag.max = 25, plot=F)
y1 <- y$acf
y2 <- c(y1)
plot(c(1,25),c(-.5,.5),xlab="",ylab="",type="n")
lines(c(1:25),y2[2:26])
abline(1.96/sqrt(l1-1),0,lty=3)
abline(-1.96/sqrt(l1-1),0,lty=3)
mtext(s=1,l=2,"lag (months)")
mtext(s=2,l=2,"autocorrelation")
mtext(s=3,l=1.5, "Figure 4:  Correlogram of first differences")
##dev.off()


### STACKING

##postscript("5to6.ps")
par(mfrow=c(1,2))

nyears <- l1/12  ### 17 years of data

plot(c(1,12),c(0,25), xlab="month", ylab="year",type="n", 
main="Figure 5: Stacked time-plot")
for(i in c(1:nyears)){
lines(c(1:12),i+x1[(12*i-11):(12*i)],lty=2)
}


plot(c(1,12),c(0,18), xlab="month", ylab="year",type="n", 
main="Figure 6: Year-adjusted stacked time-plot")
for(i in c(1:nyears)){
lines(c(1:12),i+.7*(x1[(12*i-11):(12*i)]-x1[12*i-11]),lty=i)
}

##dev.off()



### FIT TREND BY REGRESSION OR MA

index <- c(1:l1)/12
line1 <- lsfit(index,x1)

##postscript("7to10.ps")
par(mfrow=c(2,2))
plot(index, x1,xlab="year", ylab="yield(%)", type="l",lty=2)
mtext(s=3,l=1.5, "Figure 7:  Time plot with regression line")
abline(line1)

plot(index, line1$resid, xlab="year", ylab="regression residuals (%)",
type="l",lty=2)
mtext(s=3,l=1.5, "Figure 8: Linear regression residuals")


ma0 <- c(x1[12:1],x1,x1[l1:(l1-11)])
ma1 <- filter(ma0, rep(1/25,25))
ma2 <- ma1[13:(l1+12)]

plot(index, x1,xlab="year", ylab="yield(%)", type="l",lty=2)
mtext(s=3,l=1.5, "Figure 9:  time-plot with simple MA(2*12+1) filter")
lines(index, ma2)


plot(index, x1-ma2,xlab="year", ylab="MA filter residuals (%)", 
type="l",lty=2)
mtext(s=3,l=1.5, "Figure 10: MA(2*12+1) filter residuals")

##dev.off()



### partial acf
y <- acf(x2, type="partial", lag.max = 20, plot=F)
y1 <- y$acf
y2 <- c(y1)
##postscript("11to12.ps")
par(mfrow=c(1,2))
plot(c(1,20),c(-.8,.8),xlab="",ylab="",type="n")
lines(c(1:20),y2[1:20])
se <- 1.96/sqrt(l2)
abline(h=c(se,-1*se),lty=3)  
mtext(s=1,l=2,"lag (months)")
mtext(s=2,l=2,"partial autocorrelation")
mtext(s=3,l=1.5, "Figure 11:  Partial ACF of first differences")
#dev.off()


### Fit MA(1)
y <- arima0(x2,order=c(1,0,0))
y$coef
y$sigma2
y$aic
res1 <- y$resid
plot(c(0,l2/12),c(-1,1),xlab="", ylab="", type="n")
points(c(1:l2)/12,res1)
mtext(s=1,l=2,"year")
mtext(s=2,l=2,"Residuals (%)")
mtext(s=3,l=1.5, "Figure 12:  Time-plot of ARIMA (1,0,0) Residuals")

## Correlogram of MA(1) residuals
##postscript("13to14.ps")
par(mfrow=c(1,2))
y <- acf(res1, type="correlation", lag.max = 25, plot=F)
y1 <- y$acf
y2 <- c(y1)
plot(c(1,25),c(-.5,.5),xlab="",ylab="",type="n")
lines(c(1:25),y2[2:26])
abline(1.96/sqrt(l1-1),0,lty=3)
abline(-1.96/sqrt(l1-1),0,lty=3)
mtext(s=1,l=2,"lag (months)")
mtext(s=2,l=2,"autocorrelation")
mtext(s=3,l=1.5, "Figure 13:  Correlogram of ARIMA (1,0,0) Residuals")


### Fit ARMA(1,1)
y <- arima0(x2,order=c(1,0,1))
res2 <- y$resid
plot(c(0,l2/12),c(-1,1),xlab="", ylab="", type="n")
points(c(1:l2)/12,res2)
mtext(s=1,l=2,"year")
mtext(s=2,l=2,"Residuals (%)")
mtext(s=3,l=1.5, "Figure 14:  Time-plot of ARIMA (1,0,1) Residuals")
#dev.off()