# Plot some time series
par ( mfrow=c(3,2) )
ts.plot ( bonds.yield[,1], main="daily bond yields", xlab="day",
ylab="yield" )
ts.plot( co2, main="Mauna Loa CO2 concentration", xlab="Year",
ylab="concentration" )
ts.plot ( tsmatrix ( corn.rain, corn.yield ),
main="Rainfall and Corn Yield", xlab="Year",
ylab="rain (inches), yield (bushels per acre)" )
ts.plot ( lynx, main="Annual Number of Lynx Trappings", xlab="Year",
ylab="Number of Lynx" )
ts.plot ( ozone.median, main="Ozone levels",
sub="Median daily maxima for June to August", xlab="Day" )
ts.plot ( sunspots, main="Monthly Mean Relative Sunspot Numbers",
xlab="Year" )
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Fit linear and quadratic trends to the CO2 data----------------------
par ( mfrow=c(5,1) )
co2.frame _ data.frame ( as.vector(co2), as.vector(time(co2)) )
names ( co2.frame ) _ c ( "co2", "time" )
lin.fit _ lm ( co2 ~ time, co2.frame )
quad.fit _ update ( lin.fit, .~. + time^2 )
lin.pred _ predict.lm ( lin.fit, co2.frame )
quad.pred _ predict.lm ( quad.fit, co2.frame )
plot ( time(co2), co2, type="l", lty=2, xlab="Year", ylab="Concentration",
main="Mauna Loa CO2 concentration with linear and quadratic fits" )
matlines ( time(co2), cbind ( lin.pred, quad.pred ), lty=1 )
plot ( time(co2), quad.fit$res, type="l", xlab="Year", ylab="residuals" )
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Apply Moving Average filter to CO2 data------------------------------
filt _ c ( 1/24, rep(1/12,11), 1/24 )
ma.filt _ filter ( co2, filt )
tsplot ( co2, ma.filt, lty=1, xlab="Year",
main="Mauna Loa CO2 concentration with MA filter" )
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Difference the CO2 data----------------------------------------------
first.diff _ diff ( co2, differences=1 )
second.diff _ diff ( co2, differences=2 )
tsplot ( first.diff, main="First Difference of CO2 Series", xlab="Year" )
tsplot ( second.diff, main="Second Difference of CO2 Series", xlab="Year" )
#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Fit Seasonality to the CO2 Data--------------------------------------
par ( mfrow=c(4,1) )
resids _ co2 - ma.filt
resids _ matrix ( resids[!is.na(resids)], nrow=12 )
W _ apply ( resids, 1, mean )
s.hat _ W - mean(W)
plot(s.hat)
# Note that ma.filt[1:6] were NA and therefore s.hat[1] actually
# applies to observations 7, 19, 31, .... So let's rearrange.
s.hat _ s.hat[ c ( 7:12, 1:6 ) ]
# Now make s.hat into a vector the same length as co2
s.hat _ rep ( s.hat, length=length(co2) )
d _ co2 - s.hat
# Plot the deseasonalized CO2 data. A trend remains
tsplot ( d, lty=1, xlab="Year",
main="Deseasonalized Mauna Loa CO2 concentration " )
# Now reestimate the trend with MA filter. No need for a filter
# that accounts for seasonality. Or, could fit a quadratic trend
# instead of MA
filt2 _ c ( 1, 3, 5, 3, 1 )
filt2 _ filt2 / sum(filt2)
ma.filt2 _ filter ( d, filt2 )
tsplot ( co2, ma.filt2, lty=1, xlab="Year",
main="Mauna Loa CO2 concentration with MA filter
after deseasonalization" )
# Let's also plot the deseasonalized, detrended resids
tsplot ( co2 - s.hat - ma.filt2, xlab="Year",
main="Deseasonalized Mauna Loa CO2 residuals " )
#----------------------------------------------------------------------
#----------------------------------------------------------------------
#Seasonal Difference the CO2 Data--------------------------------------
par ( mfrow=c(3,1) )
season.diff _ diff ( co2, lag=12, differences=1 )
tsplot ( season.diff, main="Seasonally Differenced CO2 Series",
xlab="Year" )
# Apparently a trend remains, and perhaps a four year cycle
# Let's fit the trend with a MA filter
ma.filt2 _ filter ( season.diff, filt2 )
tsplot ( season.diff, ma.filt2, lty=c(2,1), xlab="Year",
main="Seasonally differenced CO2 and MA trend" )
# This trend is too jagged; let's make it smoother
filt3 _ rep ( 1/9, 9 )
ma.filt3 _ filter ( season.diff, filt3 )
tsplot ( season.diff, ma.filt3, lty=c(2,1), xlab="Year",
main="Seasonally differenced CO2 and MA trend" )
#-------------------------------------------------------------------
#####################################################################
par ( mfrow=c(2,2) )
wn _ rts ( rnorm(201) )
ts.plot ( wn, main="White Noise" )
acf(wn) # plot the sample autocorrelation function
# Generate an MA(1) series, estimate its autocorrelation, compare to theory
filt _ c ( 1, .5 )
MA _ filter ( wn, filt, sides=1 )
ts.plot ( MA, main="Moving Average (1)" )
acf(MA) # plot the sample autocorrelation function
#####################################################################
#####################################################################
# Generate an AR(1) series, estimate its autocorrelation, compare to theory
par ( mfrow=c(3,2) )
AR _ rep ( NA, 200 )
AR[1] _ wn[1]
for ( i in 2:200 ) AR[i] _ .5*AR[i-1] + wn[i]
ts.plot ( as.rts(AR), main="Autoregressive (1)" )
acf(AR)
#####################################################################
#####################################################################
# Let's also try some series with larger coefficients
filt _ c ( 1, .8 )
MA _ filter ( wn, filt, sides=1 )
ts.plot ( MA, main="Moving Average (1)" )
acf(MA) # plot the sample autocorrelation function
AR _ rep ( NA, 200 )
AR[1] _ wn[1]
for ( i in 2:200 ) AR[i] _ .8*AR[i-1] + wn[i]
ts.plot ( as.rts(AR), main="Autoregressive (1)" )
acf(AR)
#####################################################################
#####################################################################
# Plot the sample autocorrelation functions of our six time series
par ( mfrow=c(6,2) )
ts.plot ( bonds.yield[,1], main="daily bond yields", xlab="day",
ylab="yield" )
acf ( bonds.yield[,1] )
ts.plot( co2, main="Mauna Loa CO2 concentration", xlab="Year",
ylab="concentration" )
acf ( co2, lag.max=36 )
ts.plot ( corn.rain, main="Rainfall", xlab="Year", ylab="Inches" )
acf ( corn.rain )
ts.plot ( lynx, main="Annual Number of Lynx Trappings", xlab="Year",
ylab="Number of Lynx" )
acf ( lynx )
ts.plot ( ozone.median, main="Ozone levels",
sub="Median daily maxima for June to August", xlab="Day" )
acf ( ozone.median, lag.max=24 )
ts.plot ( sunspots, main="Monthly Mean Relative Sunspot Numbers",
xlab="Year" )
acf ( sunspots )
#####################################################################
#####################################################################
# Let's look again at the CO2 series, detrended, deseasonalized
par ( mfrow=c(2,2) )
filt2 _ c ( 1, 3, 5, 3, 1 )
filt2 _ filt2 / sum(filt2)
ma.filt2 _ filter ( d, filt2 )
# plot the original series and the estimate of trend
tsplot ( co2, ma.filt2, lty=1, xlab="Year",
main="Mauna Loa CO2 concentration with MA filter after deseasonalization" )
# plot the deseasonalized, detrended resids
tsplot ( co2 - s.hat - ma.filt2, xlab="Year",
main="Deseasonalized Mauna Loa CO2 residuals " )
# plot the acf of the residuals
acf ( co2 - s.hat - ma.filt2 )
# Tentatively, the resids look like MA(1) with a negative coefficient
#####################################################################