%pylab inline
%load_ext rmagic

%R library(astsa)
%R library(forecast)

%%R
## find linear trend in temprature
data(gtemp)
fit <- lm(gtemp ~ time(gtemp))
plot(gtemp, type="o", ylab="Global Temp Deviation")
abline(fit)
print(summary(fit))

%%R -w 500 -h 500 -u px
## Polution, Temperature and Mortality
data(cmort); data(tempr); data(part)
par(mfrow=c(3, 1))
plot(cmort, main = "Cardiovascular Motality")
plot(tempr, main = "Temperature")
plot(part, main = "particulates")


pairs(data.frame(Mortality=cmort, 
                 Temperature=tempr, 
                 Particulates=part))

fit <- lm(cmort ~ time(cmort) + tempr + tempr^2 + part)
print(summary(fit))
plot(residuals(fit), type="l")
Acf(residuals(fit))

%%R
## Regression with Lagged Variables

data(soi); data(rec)

## ccf reveals that SOI is 6 months leading of REcuritment
par(mfrow=c(1, 1))
r <- ccf(soi, rec, lag.max = 50)
## use ts.intersect to align data prior to regression
rec.soi <- ts.intersect(rec, soi.lag6 = lag(soi, -6), dframe = T)
fit <- lm(rec ~ soi.lag6, data = rec.soi)
print(summary(fit))

plot.ts(rec[-(1:6)], col = 1, lty = 1)
lines(predict(fit, rec.soi), col = 2, lty = 2)
legend("topleft", c("real", "predict"), col=1:2, lty=1:2)

%%R
## examples of nonstartionary series from real world example
data(jj)
plot(jj, xlab = "")
mtext(side = 3, line = 1, "JJ series (NONSTATIONARY): \n
Exponentially Increasing Trend and Increasing Variance")

%%R
## detrend the global temperature

## find the linear trend
fit <- lm(gtemp ~ time(gtemp))

par(mfrow = c(3, 1))
plot(gtemp, type = 'o', main = 'original gtemp')
plot(residuals(fit), type = 'o', main = 'detrended')
plot(diff(gtemp), type = 'o', main = 'first difference')

par(mfrow = c(3, 1))
Acf(gtemp, lag.max=50, main = "gtemp")
Acf(resid(fit), 50, main = 'detrended')
Acf(diff(gtemp), 50, main = "first difference")

%%R
## Acf and Scatterplot matrix of SOI
Acf(soi, lag.max=50, main = "acf finds linearly predictability only")
lag.plot(soi, 12, diag.col = "red")

%%R
## ccf and Scatterplot matrix of Rec ~ SOI_lags

%%R
## Decomposition by triangles to recover signal and noise
set.seed(1)
## assume omega(frequency) = 1/500 is known, but phase phi
## altitude A and noise level are unknown
x <- 2*cos(2*pi*(1:500)/50 + .6*pi) + rnorm(500, 0, 5)
z1 <- cos(2*pi*(1:500)/50); z2 <- sin(2*pi*(1:500)/50)
fit <- lm(x ~ z1 + z2 + 0) # zero to exclude the intercept
plot.ts(x, lty = 2)
lines(fitted(fit), lwd=2)

%%R
## use peridogram to discover signal in noise
## this time we dont know the frequency omega
set.seed(1)
x <- 2*cos(2*pi*(1:500)/50 + .6*pi) + rnorm(500, 0, 5)
I = abs(fft(x))^2/500 # the periodogram
P = (1/500)*I[1:250] # the scaled periodogram
f = (0:249)/500
plot(f, P, type = 'h')

%%R
## moving average smoother to find trend and seansonal components
x <- cmort ## weekly data
## monthly
x5 <- filter(x, sides = 2, 
             method = "convolution", filter = rep(1, 5)/5)
## yearly
x53 <- filter(x, sides = 2,
             method = "convolution", filter = rep(1, 53)/53)
plot(x, col = 1, type = 'p')
lines(x5, col = 2)
lines(x53, col = 3)