install.packages("astsa")

library("astsa")

# Fig. 1.1. Johnson & Johnson quarterly earnings per share, 84 quarters, 1960-I to 1980-IV.
plot(jj, type="o", ylab="Quarterly Earnings per Share")



# Fig. 1.2. Yearly average global temperature deviations (1880{2009) in degrees centigrade.
plot(globtemp, type="o", ylab="Global Temperature Deviations")



# Fig. 1.3. Speech recording of the syllable aaa 
 
 
 hhh sampled at 10,000 points 
# per second with n = 1020 points.
plot(speech)  



# Fig. 1.4. Returns of the NYSE. The data are daily value weighted market returns from 
# February 2, 1984 to December 31, 1991 (2000 trading days). The crash of October 19, 
# 1987 occurs at t = 938.
# djia is in astsa version 1.5+ - this is where it came from
# library(TTR)                            # install it first
# djia  = getYahooData("^DJI", start=20060420, end=20160420, freq="daily")
install.packages("xts")
library(xts)                              # if don't use TTR
djiar = diff(log(djia$Close))[-1]         # approximate returns
plot(djiar, main="DJIA Returns", type="n")  
lines(djiar) 



# Fig. 1.5. Monthly SOI and Recruitment (estimated new fish), 1950-1987.
par(mfrow = c(2,1))  # set up the graphics
plot(soi, ylab="", xlab="", main="Southern Oscillation Index")
plot(rec, ylab="", xlab="", main="Recruitment") 



# Fig. 1.6. fMRI data from various locations in the cortex, thalamus, and cerebellum; 
# n = 128 points, one observation taken every 2 seconds.
par(mfrow=c(2,1), mar=c(3,2,1,0)+.5, mgp=c(1.6,.6,0))  
ts.plot(fmri1[,2:5], col=1:4, ylab="BOLD", xlab="", main="Cortex")
ts.plot(fmri1[,6:9], col=1:4, ylab="BOLD", xlab="", main="Thalamus & Cerebellum")
mtext("Time (1 pt = 2 sec)", side=1, line=2) 



# Fig. 1.7. Arrival phases from an earthquake (top) and explosion (bottom) at 40 points per second.
par(mfrow=c(2,1))
plot(EQ5, main="Earthquake")
plot(EXP6, main="Explosion")



# Fig. 1.9. Autoregressive series generated from model (1.2).
w = rnorm(550,0,1)  # 50 extra to avoid startup problems
x = filter(w, filter=c(1,-.9), method="recursive")[-(1:50)]
plot.ts(x, main="autoregression")




# Fig. 1.10. Random walk, sigma_w = 1, with drift  delta=0.2 (upper jagged line), 
# without drift, delta=0 (lower jagged line), and a straight line with slope 0.2 (dashed line).
set.seed(154) # so you can reproduce the results
w = rnorm(200); x = cumsum(w) # two commands in one line
wd = w +.2;    xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55), main="random walk", ylab='')
lines(x, col=4) 
abline(h=0, col=4, lty=2)
abline(a=0, b=.2, lty=2)




# Fig. 1.11. Cosine wave with period 50 points (top panel) compared with the cosine wave
# contaminated with additive white Gaussian noise, sigma_w = 1 (middle panel) and 
# sigma_w = 5 (bottom panel); see (1.5).
cs = 2*cos(2*pi*(1:500)/50 + .6*pi)
w = rnorm(500,0,1)
par(mfrow=c(3,1), mar=c(3,2,2,1), cex.main=1.5)   # help(par) for info
plot.ts(cs, main = expression(x[t]==2*cos(2*pi*t/50+.6*pi)))
plot.ts(cs + w, main = expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,1)))
plot.ts(cs + 5*w, main = expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,25)))



# Fig. 1.13. ACF of the speech series.
par(mfrow=c(1,1))
acf(speech, 250)

# Fig. 1.14. Sample ACFs of the SOI series (top) and of the Recruitment series (middle), 
# and the sample CCF of the two series (bottom); negative lags indicate SOI leads 
# Recruitment. The lag axes are in terms of seasons (12 months).
par(mfrow=c(3,1))
acf(soi, 48, main="Southern Oscillation Index")
acf(rec, 48, main="Recruitment")
ccf(soi, rec, 48, main="SOI vs Recruitment", ylab="CCF")