# Script for Class 1
# Note that SS data is at www.stat.pitt.edu/stoffer/tsa2


# set working directory for convenience

setwd("/Users/bob/work/courses/stat910/")



### J&J Quarterly Earnings per Share ###

# read SS data on J&J quarterly earnings
jj.raw <- scan("data/SS_1/jj.dat")
plot(jj.raw)  
jj.ts <- ts(jj.raw, start=c(1960,1), frequency = 4)  # time labels, but come at cost
plot(jj.ts)
plot(log(jj.ts))

# Want to fit a regression on time, but regression complicates ts objects
# Simply work with regular R series

jj      <- jj.raw
time    <- 1960+(0:83)/4
quarter <- as.factor(rep(1:4,21))  # quarter labels

plot(time,log(jj), main="Log of J&J Quarterly Earnings Per Share")

# Fit a regression on time, residuals have patterns
regr <- lm(log(jj) ~ time); summary(regr)
abline(regr, col="red")   # add to plot

plot(time, residuals(regr), col=1:4)  # color is nice for this

# Revise regression by adding quarterly indicators
regr.2 <- lm(log(jj) ~ time + quarter); summary(regr.2)
plot(time, residuals(regr.2), col=1:4)  # color is nice for this

regr.3 <- lm(log(jj) ~ time + quarter + quarter*time); summary(regr.3)
plot(time, residuals(regr.3), col=1:4)  # color is nice for this

# Q: What's happening?  Certainly not random pattern left in residuals.
# Q: What's the presence of a pattern here say about the regressions?


### Global temperature deviations ###

temp <- scan("data/nasa_global.dat")
time <- 1880 + (0:(length(temp)-1))/12
plot(time,temp, type="l", main="NASA Global Temperature Deviation",
     xlab="Year", ylab="1/100 Degree Deviation")

# plot deviations of deviations
smth <- loess(temp ~ time)
lines(time,smth$fitted, col="red")

smth.4 <- loess(temp ~ time, span=0.4)
lines(time,smth.4$fitted, col="blue")

plot(time, smth.4$residuals)



### Speech Signal ###

# Much longer, with quite different and regular repeating pattern

speech <- scan("data/SS_1/speech.dat")
length(speech)
plot(speech, type="l")

# look at the autocorrelations
acf(speech)

# look at the periodogram
spec.pgram(speech, demean=TRUE)

s <- spec.pgram(speech, taper=0.5, demean=TRUE, plot=FALSE, spans=c(25,13))
lines(s$freq, s$spec, col="red")



### Southern Oscillation and Fish ###

# Arranged as two series, with 453 months in each

soi  <- scan("data/SS_1/soi.dat")
fish <- scan("data/SS_1/recruit.dat")

par(mfrow=c(2,1))  
  plot(soi, main="Southern Oscillation", type="l")
  plot(fish, main="Fish Recruitment", type="l")
par(mfrow=c(1,1))

# no contemporaneous correlation
plot(soi,fish)

# dependence at other lags?
plot(soi[-453], fish[-1])
plot(soi[-(453-(0:3))], fish[-(1:4)])
plot(soi[-(453-(0:7))], fish[-(1:8)])

# this function computes the correlations at varying lags
ccf(soi,fish)


### Brain Imaging ###

mri <- read.table("data/SS_1/fmri.dat")
time <- 2*(1:128)

# these are very much in phase
 plot(time,mri[,2],type="l",col="red")
lines(time,mri[,3],type="l",col="blue")
lines(time,mri[,4],type="l",col="green")
lines(time,mri[,5],type="l",col="black")
lines(time,mri[,6],type="l",col="orange")

# these are rather different
 plot(time,mri[,7],type="l",col="red")
lines(time,mri[,8],type="l",col="blue")
lines(time,mri[,9],type="l",col="green")



### Earthquake or Explosion ###

# Arranged as 2048 of earthquake, then explosion (sampled at 40/second)
boom <- matrix(scan("data/SS_1/eq5exp6.dat"), ncol=2)
length(boom[,1])
par(mfrow=c(2,1))  
  plot(boom[,1], main="Earthquake", type="l")
  plot(boom[,2], main="Explosion", type="l")
par(mfrow=c(1,1))



### Random Walk ###

# can be very hard to distinguish this sort of noise from results of a model
# since you only have one realization, one sample path
rw <- ts(cumsum(rnorm(100))); plot(rw)


### Linear Process ###
xt <- filter(rnorm(100), filter=c(.2,.4,.4,.4,.2), method="convolution"); plot(xt)


### Autoregression ###
xt <- filter(rnorm(100), filter=c(.2,.5,.2), method="recursive"); plot(xt)


# How would you distinguish among these?