# Plot some timeseries
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" )
#----------------------------------------------------------------------
# Work with co2 first
#
# Apply Moving Average filter to CO2 data------------------------------
par ( mfrow=c(4,1) )
filt _ c ( 1/24, rep(1/12,11), 1/24 )
ma.filt _ filter ( co2, filt )
ts.plot ( co2, ma.filt, lty=1, xlab="Year",
main="Mauna Loa CO2 concentration with MA filter" )
resids _ co2 - ma.filt
ts.plot ( resids, lty=1, xlab="Year",
main="Mauna Loa CO2 residuals from MA filter" )
filt2 _ c ( 1, 3, 5, 3, 1 )
filt2 _ filt2 / sum(filt2)
ma.filt2 _ filter ( d, filt2 )
ts.plot ( co2, ma.filt2, lty=1, xlab="Year",
main="Mauna Loa CO2 concentration with MA filter after deseasonalization")
resids2 _ co2 - ma.filt2
ts.plot ( resids2, lty=1, xlab="Year",
main="Mauna Loa CO2 residuals from MA filter" )
par ( mfrow=c(3,1) )
# Estimated spectrum from data with trend. Note the concentration at low
# frequency
spectrum ( co2, detrend=F ) # without smoothing
p2 <- spec.pgram( co2, spans=c(9,7), detrend=F, plot=F ) # with smoothing
spec.plot(p2,add=T)
# Estimated spectrum from smoothed data
spectrum ( ma.filt, detrend=F ) # without smoothing
p2 <- spec.pgram( ma.filt, spans=c(9,7), detrend=F, plot=F ) # with smoothing
spec.plot(p2,add=T)
spectrum ( ma.filt2, detrend=F ) # without smoothing
p2 <- spec.pgram( ma.filt2, spans=c(9,7), detrend=F, plot=F ) # with smoothing
spec.plot(p2,add=T)
# Estimated spectrum of detrended data
spectrum ( co2 ) # without smoothing
p2 <- spec.pgram( co2, spans=c(9,7), plot=F ) # with smoothing
spec.plot(p2,add=T)
# Estimated spectrum from smoothed data
spectrum ( ma.filt ) # without smoothing
p2 <- spec.pgram( ma.filt, spans=c(9,7), plot=F ) # with smoothing
spec.plot(p2,add=T)
spectrum ( ma.filt2 ) # without smoothing
p2 <- spec.pgram( ma.filt2, spans=c(9,7), plot=F ) # with smoothing
spec.plot(p2,add=T)
# Estimated spectrum of residuals
# Estimated spectrum from smoothed data
spectrum ( resids ) # without smoothing
p2 <- spec.pgram( resids, spans=c(9,7), plot=F ) # with smoothing
spec.plot(p2,add=T)
spectrum ( resids2 ) # without smoothing
p2 <- spec.pgram( resids2, spans=c(9,7), plot=F ) # with smoothing
spec.plot(p2,add=T)
# transfer function
f <- p2$freq
f <- f[f
tf <- rep(NA,length(f))
tf2 <- rep(NA,length(f))
tf3 <- rep(NA,length(f))
for ( i in seq(along=f) ) {
tf[i] _ 1 - 1/12 - 2 * sum ( filt[1:6]*cos((1:6)*f[i]) )
tf2[i] _ 1 - 5/13 - 2 * ( (1/13)*cos(f[i]) + (3/13)*cos(2*f[i]) )
tf3[i] _ 1 - 2*(1/3)*cos(f[i])
}
matplot(f,cbind(tf^2,tf2^2,tf3^2),type="l",lty=1:3)
#----------------------------------------------------------------------