# TIME SERIES ANALYSIS
# https://www.statology.org/ / good source

# set working directory and clear memory
setwd("C:/data/R_data/GGEC")
rm(list = ls())

# read from the text file 
bt <- read.table("Brno_airport_mean_air_temperature.csv", header = T, dec = '.', sep = ';')


###############################
# bt <- read.table("Ireland_annual_temerature.csv", header = T, dec = '.', sep = ';')
# bt <- read.table("Slovenia_annual_temerature.csv", header = T, dec = '.', sep = ';')
# bt <- read.table("Zambia_annual_temerature.csv", header = T, dec = '.', sep = ';')
# Downloaded from Berkeley earth: https://berkeleyearth.org/data/
# Data are anomalies w.r.t 1951 - 1980
###############################

names(bt)
head(bt)
tail(bt)
dim(bt)


# Always start with EDA (Exploratory Data Analysis)
# simple plot
plot(bt$YEAR, bt$ANN)
plot(bt$YEAR, bt$ANN, type="l")
summary(bt$ANN)
hist(bt$ANN)

# qq-plot
qqnorm(bt[,2])
qqline(bt[,2])

# normality test
shapiro.test(bt$ANN)

# create time series object
temperature <- ts(bt$ANN, start=bt$YEAR[1])
str(temperature) # shows object structure
plot.ts(temperature, lw=2, col="blue")

# find maximum teperature and year of ocurrence
temperature[which.max(temperature)]
time(temperature)[which.max(temperature)]


# Subset
bt61 <- window(temperature, start = 1961,end = 1990)
bt91 <- window(temperature, start = 1991,end = 2019)

# set space for two graphs
par(mfrow=c(1,2),las=1)
plot.ts(temperature,lty=4,,lw=1, col="black")
lines(bt61,lty=1,,lw=2, col="blue")
lines(bt91,lty=1,,lw=2, col="red")

boxplot(bt61, bt91, col=c("blue", "red"))
t.test(bt61, bt91)

################################################################################################
# TREND ANALYSIS / 1. Linear regression
# parameters estimated with the Least Squares Method

# Simple Linear regression Model trend.1 <- lm(bt$ANN ~ bt$YEAR) 
# another possibility  trend.1 <- lm(temperature ~ time(temperature))
# To save repeteated writings we use:
x <- as.vector(bt$YEAR)
y <- as.vector(bt$ANN)
trend.1 <- lm(y ~ x)
summary(trend.1)

# plot(temperature)
par(mfrow=c(1,1))
plot(bt$YEAR, bt$ANN, type="l",xlab="Year", ylab="Temperature anomaly")
abline(trend.1, lw=3, col="blue")

# EXPLORE STRUCTURE OF RESULTS - very useful to find what You want
str(trend.1)
trend.1$coefficients
slope <- trend.1$coefficients[2]
slope

# TREND ANALYSIS / 2. NONPARAMETRIC TREND ESTIMATE 
# 1) Mann-Kendal test for monotonic trend 
# 2) Theil-Sen parameters estimate 
# https://rcompanion.org/handbook/F_12.html

# we nedd special packages
if(!require(Kendall)){install.packages("Kendall")}
library(Kendall)

# Mann-Kendal test
MKD <- MannKendall(y)
# Interpretation using p-value
summary(MKD)

if(!require(mblm)){install.packages("mblm")}
library(mblm)

trend.2 = mblm(y ~ x)
summary(trend.2)
# Add confidence limits 
confint(trend.2)

abline(trend.2, lw=3, col="red")

# Add legend
legend(x="bottomright", legend=c("Least Square Regression", "Non-parametric"), bty="n", lwd=3, col=c("blue", "red"))
#####################################################################################

# TREND MODEL VERIFICATION

# plot resiuals
res <- resid(trend.1)
plot(fitted(trend.1), res)
abline(0,0)

#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 

# calculate RMSE of regression model
sqrt(mean(trend.1$residuals^2))

#####################################################################################

# ADAPTIVE TREND - 1. MOVING AVERAGES
plot(x, y, type="l",xlab="Year", ylab="Temperature", main=" Moving averages")

# Moving averages
if(!require(forecast)){install.packages("forecast")}
library(forecast)

sm.1 <- ma(y, order=5)
lines(x,sm.1, col="blue",lwd=3)
sm.1 <- ma(y, order=15)
lines(x,sm.1, col="green",lwd=3)
sm.1 <- ma(y, order=30)
lines(x,sm.1, col="red",lwd=3)

# Add legend
legend(x="bottomright", legend=c("ma_05", "ma_15", "ma_30"), bty="n", lwd=3, col=c("blue", "green", "red"))


##########################################################################################
# ADAPTIVE TREND - 2. GAUSSIAN LOW PASS FILTER

if(!require(smoother)){install.packages("smoother")}
library(smoother)

plot(x, y, type="l",xlab="Year", ylab="Temperature", main=" Gaussian filter")

sm.2 <- smth.gaussian(y, window=10, tails=TRUE)
lines(x,sm.2, col="magenta",lwd=3)

sm.2 <- smth.gaussian(y, window=20, tails=TRUE)
lines(bt$YEAR,sm.2, col="red",lwd=3)

sm.2 <- smth.gaussian(y, window=30, tails=TRUE)
lines(x,sm.2, col="cyan",lwd=3)

legend(x="bottomright", legend=c("gf_10", "gf_20", "gf_30"), bty="n", lwd=3, col=c("magenta", "red", "cyan"))
#############################################################################################
# ANOTHER APPROACHES TO SMOOTHING
# Local regression
# https://www.statology.org/lowess-smoothing-r/

sm.3 <- lowess(y  ~ x, f=1/5)
plot(x, y, type="l",xlab="Year", ylab="Temperature", main="Time series smoothing approaches")
lines(x,sm.3$y, col="darkgreen",lwd=3)


# SPLINES 
# http://users.stat.umn.edu/~helwig/notes/smooth-spline-notes.html
lines(smooth.spline(x, y, df=18),col="magenta",lwd=3)
#############################################################################################