#########################################################################
### Illustration of smoothing cardiovascular mortality series from
### Shumway & Stoffer (2006) [http://www.stat.pitt.edu/stoffer/tsa2/textRcode.htm]
### The data are average weekly cardiovascular mortality in Los Angeles
### County over the 10-year period 1970-1979.
#########################################################################
# read data and convert to a time series object
mort = ts(scan("Data/cmort.dat"),start=1970, frequency=52) 
#
# Moving Average smooothers: 
ma5 = filter(mort, sides=2, rep(1,5)/5)
ma53 = filter(mort, sides=2, rep(1,53)/53)
plot(mort, type="p", xlab="week", ylab="mortality")
lines(ma5)  # yearly smooth (trend)
lines(ma53) # monthly smooth (trend + seasonality)
#
# Polynomial and Harmonic Regression smoothers: 
wk = time(mort) - mean(time(mort))  # week (centered) - wk is basically t/52 ...
wk2 = wk^2                          # ... because frequency=52 for mort
wk3 = wk^3
c = cos(2*pi*wk)  
s = sin(2*pi*wk)
reg1 = lm(mort~wk + wk2 + wk3, na.action=NULL)
reg2 = lm(mort~wk + wk2 + wk3 + c + s, na.action=NULL)
plot(mort, type="p", xlab="week", ylab="mortality")
lines(fitted(reg1)) # trend
lines(fitted(reg2)) # trend + seasonality
#
# Kernel Smoothing:
plot(mort, type="p", xlab="week", ylab="mortality")
lines(ksmooth(time(mort), mort, "normal", bandwidth=2)) # trend
lines(ksmooth(time(mort), mort, "normal", bandwidth=10/52)) # trend + seasonality 
#
# Lowess smoothing:
plot(mort, type="p", xlab="week", ylab="mortality", main="lowess")
lines(lowess(mort, f=.02)) # trend + seasonality    
lines(lowess(mort, f=2/3)) # trend
#
# Splines smoothing:
plot(mort, type="p", xlab="week", ylab="mortality")
lines(smooth.spline(time(mort), mort, spar=0)) # trend + seasonality         
lines(smooth.spline(time(mort), mort, spar=1)) # trend