y=scan('D:/334/CumbriaTempMonthly.txt',skip=2,nlines=240)

#remove the last 10% of data and plot the time series
y=scan('D:/334/CumbriaTempMonthly.txt',skip=2,nlines=216)
yts=ts(y)
date=seq(2000,2017+11/12,1/12)
plot(date,yts,main="MONTHLY TEMP IN CUMBRIA",type='p')

#linear regression model
n=length(y)
time=seq(1:n)#^1.15
# where n=216 and period=12
month=c(rep(seq(1:12),18))
fmonth=as.factor(month)
fmonth

ymod1=lm(y~time)
summary(ymod1)
yfit1=ymod1$fitted
lines(date,yfit1,col=2)

#Residuals of linear regression - use to check independence assumption
yres1=ymod1$res
plot(date,yres1)
lines(date,yres1)
#are they independent? - lecture notes

#plotting sample autocorrelation and partial autocorrelation function of residuals
acf(yres1)
pacf(yres1)
#numerical values
u=acf(yres1)

#estimate spectral density - periodogram of residuals
yspec=spec.pgram(yres1,pad=4,detrend=TRUE)
plot(yspec$freq,yspec$spec,type='l',col=5)

# harmonic regression
FF = 1/12 # Fundamental frequency
ymod2=lm(y~time+cos(2*pi*time*FF)+ sin(2*pi*time*FF))
summary(ymod2)
yfit2= ymod2$fitted
plot(y, main="Temperature and harmonic regression (fundamental only)",type='p',col=2)
lines(yfit2,col=4)
lines(yfit1,col=3)

# residuals of harmonic regression
ts.plot(ymod2$res, main="Residuals after harmonic regression (fundamental only)",col=4)

# periodogram of residuals
yspec=spec.pgram(ymod2$res,pad=4,main="Periodogram of residuals after harmonic regression (fundamental only)")
plot(yspec$freq,yspec$spec,type='l',col=4,main="periodogram of residuals after harmonic regression (fundamental only)")

# acf and pacf of harmonic residuals
acf(ymod2$res,main="acf of residuals after harmonic regression (fundamental only)")
pacf(ymod2$res,main="acf of residuals after harmonic regression (fundamental only)")

# Ljung Box statistics
#original
p = rep(0, 15)
for (i in 1:15) {
  p[i] = Box.test(ymod1$res, i, type = "Ljung-Box" )$p.value
}
p.min = .05 #level of significance
plot(p,ylim=c(0,1),main="p values for Ljung-Box statistic", xlab="lag",ylab="p value")
abline(h=p.min,col=4,lty=2)
#harmonic
p = rep(0, 15)
for (i in 1:15) {
  p[i] = Box.test(ymod2$res, i, type = "Ljung-Box" )$p.value
}
p.min = .05 #level of significance
plot(p,ylim=c(0,1),main="p values for Ljung-Box statistic from harmonic regression", xlab="lag",ylab="p value")
abline(h=p.min,col=4,lty=2)

#harmonic regression, fundamental and 2nd harmonic
ymod3=lm(y~time+cos(2*pi*time*FF)+ sin(2*pi*time*FF)+
           cos(4*pi*time*FF)+ sin(4*pi*time*FF))
summary(ymod3)
yfit3= ymod3$fitted
plot(y, main="Temperature and harmonic regression - second harmonic)",type='p',col=2)
lines(yfit3,col=4)
#residuals of harmonic regression
ser=ts(ymod3$res)
plot(ser,main="Resdiuals after harmonic fit - 2nd harmonic)",col=4)

#periodogram of residuals
yspec=spec.pgram(ymod3$res,pad=4,main="periodogram of residuals after harmonic regression - 2nd harmonic, linear scale)"
plot(yspec$freq,yspec$spec,type='l',col=4,main="periodogram of residuals after harmonic regression - with 2nd harmonic, log scale)")
# acf and pacf of residuals
acf(ymod3$res, main="acf of residuals after harmonic regression 2nd harmonic")
pacf(ymod3$res, main="acf of residuals after harmonic regression 2nd harmonic")
# Ljung Box statistics 2nd harmonic
p = rep(0, 15)
for (i in 1:15) {
  p[i] = Box.test(ymod3$res, i, type = "Ljung-Box" )$p.value
}
p.min = .05  # level of significance
plot(p, ylim=c(0,1),main="p values with Ljung-Box for residuals from harmonic regression 2nd harmonic", xlab="lag",ylab="p value")
abline(h=p.min,col=4,lty=2)

#harmonic regression, 3rd harmonic
ymod4=lm(y~time+cos(2*pi*time*FF)+ sin(2*pi*time*FF)+
           cos(4*pi*time*FF)+ sin(4*pi*time*FF)+
           cos(6*pi*time*FF)+ sin(6*pi*time*FF))
summary(ymod4)
yfit4= ymod4$fitted
plot(y, main="Cumbria Temperature and harmonic regression (with 3rd harmonic)",type='p',col=2)
lines(yfit4,col=4)
# residuals of harmonic regression
ser=ts(ymod4$res)
plot(ser,main="Resdiuals after harmonic fit - 3rd harmonic")
# periodogram of residuals
yspec=spec.pgram(ymod4$res,pad=4,main="periodogram of residuals after harmonic regression - 3rd harmonic, linear scale")
plot(yspec$freq,yspec$spec,type='l',col=4, main="periodogram of residuals after harmonic regression - 3rd harmonic, log scale")
# acf and pacf of residuals
acf(ymod4$res, main="acf of residuals after harmonic regression (with 3rd harmonic)")
pacf(ymod4$res, main="pacf of residuals after harmonic regression (with 3rd harmonic)")
# Ljung Box statistics
p = rep(0, 13)
for (i in 1:13) {
  p[i] = Box.test(ymod4$res, i, type = "Ljung-Box" )$p.value
}
p.min = .05  # level of significance
plot(p, ylim=c(0,1),main="p values with Ljung-Box for residuals from harmonic regression 3 harmonic", xlab="lag",ylab="p value")
abline(h=p.min,col=4,lty=2)

#use periodogram to detect cycles in a time series
#peak in periodogram = cycle associated with a strong seasonal pattern
#remove trend before examining its periodogram
y=scan('D:/334/CumbriaTempMonthly.txt',skip=2,nlines=216)
n=length(y)
time1=seq(1:n)^1.15
#time1 is super linear whilst time is linear trend
time=seq(1:n)
ymod1=lm(y~time1)
yres1=ymod1$res
ser=ts(yres1)
plot(ser)