Time Series Analysis for Thickness of Ozone data
Introduction
Atmosphere of the earth is divided in to four layers in the below altitudes, namely Troposphere, Stratosphere, Mesosphere and Thermosphere (Lu,2015). Stratosphere which mainly consisted of Ozone (O3) and Oxygen (O2) absorb UV radiations (Dessler, 2000). This layer of Ozone is important to continue to maintain the ideal environment for animals on Earth. This study is focused on analyzing the changes of Ozone layer from 1927 to 2016. Thickness of Ozone layer has measured in Dobson Units. Aimed to forecast changes in thickness of Ozone layer for next five years starting from 2017 at the end of the study using a most fitted time series model. The data set is consisted of only one variable, the thickness of ozone layer given in Dobson units and data series is given by data1.csv. Analysis is done using functions in Rstudio and package used is “TSA”.
library(TSA)#load the TSA package ##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tar#Reading the data
Ozone <- read.table("~/Downloads/data1.csv", quote="\"", comment.char="",header = FALSE)
#converting the data frame in to a time series using ts function
Ozone <- ts(as.vector(Ozone), start = 1927, end = 2016, frequency = 1)
class(Ozone)## [1] "ts"Time Series Plot
plot(Ozone,type='o',ylab="Thickness of Ozone layer (Dobson Units)", xlab="Year", main="Time Series Plot for Thickness of Ozone", col="blue")
Time Series plot of the Ozone thickness series clearly visualise downward trend and it is obvious that positive Dobson values responsible for the increases in thickness while negative values for the decreases.Seasoanlity is not observed and behaviour of the series is Moving Average and there are no intervention points in the series. Changing varainces is not observed in the series.
Plotting the scatter plot
plot(y=Ozone, x=zlag(Ozone),ylab = "Thickness of the Ozone layer (Dobson Units)", xlab="Previous year Ozone Thicnkess", main = "Scatter plot of neighboring thickness", col="blue")
#Correlation
y=Ozone # Read the data in to y
x=zlag(Ozone) #Generate first lag of the series
index = 2:length(x)
cor(y[index],x[index]) #Calculate the correlation ## [1] 0.8700381Scatter plot visualise a strong positive correlation between the neighboring points where by the value of correlation which is 0.870 it is proven that strong positive relationship between previous year thickness of ozone layer.
Linear Trend Model and it’s Redisual Analysis
Linear Trend Model
Linear trend model is fitted to the data first,
#Linear trend model
model.Ozone1 = lm(Ozone~time(Ozone)) #label the model as model.Ozone1
summary(model.Ozone1)##
## Call:
## lm(formula = Ozone ~ time(Ozone))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7165 -1.6687 0.0275 1.4726 4.7940
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 213.720155 16.257158 13.15 <2e-16 ***
## time(Ozone) -0.110029 0.008245 -13.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.032 on 88 degrees of freedom
## Multiple R-squared: 0.6693, Adjusted R-squared: 0.6655
## F-statistic: 178.1 on 1 and 88 DF, p-value: < 2.2e-16Slope and the intercept of the linear trend model is -0.110029 and 213.72 respectively. Slope of the linear trend model is statistically significant at 5% level of significant.
p-value in the linear model is less than 0.05 which indicates that model is significant at the level of 5% while Adjusted R squared value explains that 66.55% of the variation in ozone thickness data series is explained by linear model.
Plot for Linear trend model
plot(Ozone,type='o',ylab='Thickness of the Ozone layer (Dobson units)', xlab="Year",main="Fitted least curve to Thickness of Ozone data",lty=4,col="blue")
abline(model.Ozone1)
Residual Analysis of Linear Trend Model
Since linear model showed a partial significance the residuals were saved and analyzed to receive a better understanding.
Plot for time Series for linear trend residuals
res.model.Ozone1 = rstudent(model.Ozone1) #Label residual model as res.model.Ozone1
plot(y=res.model.Ozone1,x=as.vector(time(Ozone)),xlab = 'Year', ylab = 'Standardized Residuals', main = "Time Series plot of residuals",type='o', col="blue", ylim = c(-3,3))
Downward trend of the series has improved after linear modelling of the series and it is clearly evident in the plot that series now shows less trend than the original time series plot.
Time series plot of Standardised Residuals
plot(res.model.Ozone1, x=as.vector(fitted(model.Ozone1)),xlab="Fitted Trend Values", ylab="Standardized Residuals",type='n',main="Time Series plot of Standardized Residuals",col="blue")
points(y=res.model.Ozone1,x=as.vector(fitted(model.Ozone1)))
Histogram for residuals of linear trend model
hist(res.model.Ozone1, xlab='Standardized Residuals',main="Histogram for residuals of Linear Treand Model", col = "light blue",xlim=c(-3,3))
Normality of the residuals can be explained using histogram. The histogram indicates slight symmetry in the ends of the graph and data are laid with in -3 to 3.
Normal qq Plot
qqnorm(res.model.Ozone1,col="blue")
qqline(res.model.Ozone1,col="red",lwd=2)
Shapiro-Wilk Normality Test
shapiro.test(res.model.Ozone1)##
## Shapiro-Wilk normality test
##
## data: res.model.Ozone1
## W = 0.98733, p-value = 0.5372Shapiro-Wilk test shows p-value of 0.5372 since null hypothesis is not rejected which state that the stochastic component of the linear model is normally distributed.
Sample Auto Correlation Function
acf(res.model.Ozone1,main="ACF of Standardized Residuals", col="red")
Correlations values higher than the confidence boundries are observed at several lags since stochastic component of the series is not white noise. There are three significant lags in the linear trend model which also prove that series is not a white noise.
Quadratic Trend Model and it’s Residual Analysis
Ozone thickness time series data are then modelled using quadratic trend model,
Quadratic Trend Model
t=time(Ozone)
t2=t^2
model.Ozone2=lm(Ozone~t+t2) #label the model as model.Ozone2
summary(model.Ozone2)##
## Call:
## lm(formula = Ozone ~ t + t2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1062 -1.2846 -0.0055 1.3379 4.2325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.733e+03 1.232e+03 -4.654 1.16e-05 ***
## t 5.924e+00 1.250e+00 4.739 8.30e-06 ***
## t2 -1.530e-03 3.170e-04 -4.827 5.87e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.815 on 87 degrees of freedom
## Multiple R-squared: 0.7391, Adjusted R-squared: 0.7331
## F-statistic: 123.3 on 2 and 87 DF, p-value: < 2.2e-16p-value of the quadratic model is less than 0.05 since we can assume that the model is significant at the level of 5% significant level.
R squared value of the model is 0.733 since 73.31% of the variation in the ozone layer thickness series is explained by the quadratic trend which indicates that quadratic model is more significant and fitting than the linear model that has only 66%.
Residual Analysis of Quadratic Trend Model
Time Series plot for quadratic trend residuals
res.model.Ozone2 = rstudent(model.Ozone2) #Label residual model as res.model.Ozone2
plot(y=res.model.Ozone2,x=as.vector(time(Ozone)),xlab = 'Year', ylab = 'Standardized Residuals', main = "Time Series plot of residuals",type='o', col="blue", ylim = c(-3,3))
Downaward trend in the orginial series has being removed after modelling the original data in to a quadratic model and improved than the linear model.
Time Series plot of Standardised Residuals
plot(res.model.Ozone2, x=as.vector(fitted(model.Ozone2)),xlab="Fitted Trend Values", ylab="=Standardized Residuals",type='n',main="Time Series plot of Standardized Residuals")
points(y=res.model.Ozone2,x=as.vector(fitted(model.Ozone2)))
Histogram for residuals of linear trend model
hist(res.model.Ozone2, xlab='Standardized Residuals',main="Histogram for residuals of Quadratic Treand Model", col = "light blue",xlim=c(-3,3))
Normal qq Plot
qqnorm(res.model.Ozone2)
qqline(res.model.Ozone2,col="red",lwd=2)
Shapiro-Wilk Normality Test
shapiro.test(res.model.Ozone2)##
## Shapiro-Wilk normality test
##
## data: res.model.Ozone2
## W = 0.98889, p-value = 0.6493Shapiro-Wilk test shows p-value of 0.6493 since null hypothesis is not rejected which state that the stochastic component of the quadratic model is normally distributed.
Sample Auto Correlation Function
acf(res.model.Ozone2,main="ACF of Standardized Residuals", col="red")
Four significant lags are observed after quadratic modelling and correlations values are higher than the confidence bounds at observed several lags since expected smoothness of the time series plot is not found in the white noise process.
Harmonic Trend Model and it’s Residual Analaysis
Harmonic Trend Model
Harmonic model
har.=harmonic(Ozone,0.45)
model.Ozone3=lm(Ozone~har.)
summary(model.Ozone3)##
## Call:
## lm(formula = Ozone ~ har.)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3520 -1.8906 0.4837 2.3643 6.4248
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.970e+00 4.790e-01 -6.199 1.79e-08 ***
## har.cos(2*pi*t) NA NA NA NA
## har.sin(2*pi*t) 5.462e+11 7.105e+11 0.769 0.444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.522 on 88 degrees of freedom
## Multiple R-squared: 0.006672, Adjusted R-squared: -0.004616
## F-statistic: 0.5911 on 1 and 88 DF, p-value: 0.4441Plot for Fitted Harmonic Trend Model
plot(ts(fitted(model.Ozone3)), ylim = c(min(c(fitted(model.Ozone3),as.vector(Ozone))), max(c(fitted(model.Ozone3),as.vector(Ozone)))),ylab="Thickness of Ozone layer (Dobson Units)",xlab="Year", main = "Fitted harmonic trend model to Thickness of Ozone data", type="l",lty=2,col="red")
lines(as.vector(Ozone),type="o")
Seasonality is not observed in the data and the p-value of the model is more than 0.05 suggests that model is insignificant. at the level of 5% significant level. Adjusted R squared value is also less than the values of linear and quadratic models. Missing values in the model serve as another reason for insignificance of the model.
Since the Harmonic model is insignificant a residual analysis for the model will not be conducted.
Forecasting
Carefully considering and analyzing the three models, can suggest that the quadratic model is closest with the time series of actual data hence we used quadratic model to forecast the changes in the thickness of ozone layer for next five years (2017-2021).
t = c(2017, 2018, 2019, 2020, 2021)
t2 = t^2
new = data.frame(t,t2)
forecast = predict(model.Ozone2, new, interval = "prediction")
print(forecast)## fit lwr upr
## 1 -10.34387 -14.13556 -6.552180
## 2 -10.59469 -14.40282 -6.786548
## 3 -10.84856 -14.67434 -7.022786
## 4 -11.10550 -14.95015 -7.260851
## 5 -11.36550 -15.23030 -7.500701plot(Ozone, xlim = c(1925,2023), ylim = c(-15, 5), type="o", ylab="Thickness of Ozone Layer in (Dobson Units)" , main = " Forceasting using quadratic trend model", sub="Thickness of Ozone layer from 2017 to 2021",col="blue")
lines(ts(as.vector(forecast [,1]), start = c(2017,1), frequency = 1), col="red", type="l",lwd=2)
lines(ts(as.vector(forecast [,2]), start = c(2017,1), frequency = 1), col="darkgreen", type="l",lwd=2)
lines(ts(as.vector(forecast [,3]), start = c(2017,1), frequency = 1), col="darkgreen", type="l",lwd=2)
Conclusion
Thickness of the ozone layer data are showing only downward trend while it shows no seasonality, changing variances or intervention points. After the descriptive analysis of the data, linear, quadratic and harmonic modeling are performed for the data.
According to the analysis it is evident that quadratic trend model follows the downward trend of original time series data than linear and harmonic models. Forecasting for the ozone layer data are performed using quadratic trend model.
References
Dessler, A., 2000. Chemistry and Physics of Stratospheric Ozone (Vol. 74). Elsevier.
Qing-bin, L., 2015. New theories and predictions on the ozone hole and climate change. World Scientific.pp.07