Problem

The purpose of this assignment is to analyse the yearly changes in the thickness of Ozone layer from 1927 to 2016 in Dobson units. A negative value in the dataset represents a decrease in the thickness and a positive value represents an increase in the thickness. The task is to find the good fitting trend model to the given dataset and provide a 5-step ahead forecast.

Data Description

The data is in csv format and represents 90 observations from 1927 to 2016

Descriptive Analysis

par(mar=c(5,4,4,2),cex.main=1, cex.lab=1, cex.axis=1)
#plot(ozoneLayer, type='o', ylab='Ozone layer Thickness')
plot(ozoneLayer,ylab='Thickness of Ozone layer (Dobson units)',xlab='Year',type='o', main = "Time series
plot of changes in Ozone Layer Thichness.")

Time Series Plot

Scatter Plot

From the above time series plot(Fig 1.) & scatter plot(Fig 2.) we can discuss on the following fundamental concepts mentioned below:

Trend: Clear evidence of negative trend, signifying that the thickness of ozone layer has reduced over given period of time.
Seasonality: No repeating pattern seen, signifying no apparent indication of seasonality effect
Changing variance: slight effect of changing variance seen at early 1950s and large variance in early 2000s
Autoregressive behaviour or Moving Average: the succeeding points are related to one another as the data points are affecting the next neighbouring years.
Intervention: No apparent Intervention point noticed

y = ozoneLayer
x = zlag(ozoneLayer)

excludeNAs = 2:length(x)

PearsonCor = cor(y[excludeNAs], x[excludeNAs])

knitr::kable(PearsonCor, caption = "Correlation") %>%
  column_spec(1, bold = T, color = "white", background = "#FF4040")

Correlation
x
0.8700381

As we observed from the scatter plot, the correlation factor is also supporting our evidence that there is high correlation between succeding data points. i.e. the ozone layer thickness of one year impacts the next years ozone layer thickness.

Model

Linear Regression Model

## 
## Call:
## lm(formula = ozoneLayer ~ time(ozoneLayer))
## 
## 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(ozoneLayer)  -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-16

plot(ozoneLayer,type='o',ylab='y',main = "Time series plot for simulated random walk series")
abline(lmModel, col = 'red', lty=1)

Fitted Linear Curve to Ozone Layer Thickness

The linear regression model results signifies that the intercept and the data points are significant. Also the model is also significant with p value < 0.05 with Residual standard error is also small. There is not much difference between the R². & the Adjusted R². values. From Fig 3. we can also see that the model is not a that good a fit for the given data points. Hence let us try exploring quadratic models:-

Quadratic Models

t= time(ozoneLayer)
t2 = t^2
quadModel_2 = lm(ozoneLayer~t+t2)

## 
## Call:
## lm(formula = ozoneLayer ~ 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-16

plot(ts(fitted(quadModel_2)), 
       ylim = c(min(c(fitted(quadModel_2), as.vector(ozoneLayer))),
                max(c(fitted(quadModel_2), as.vector(ozoneLayer)))),
       ylab='Ozone Layer Thickness in Dobson units', main = "Fitted quadratic curve to the data", col=c('red'))
 
lines(as.vector(ozoneLayer),type="o", col=c("black"))
legend("bottomleft",lty=1, bty = "n" ,text.width = 15, col=c("black","red"), 
       c("Ozone Layer Thickness", "Fitted Quadratic Model of Order 2"))

Fitted Quadratic Curve to Ozone Layer Thickness of Order 2

t3 = t^3
quadModel_3 = lm(ozoneLayer~t+t2+t3)

## 
## Call:
## lm(formula = ozoneLayer ~ t + t2 + t3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3466 -1.3052  0.0741  1.3259  4.4073 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.073e+05  1.064e+05   1.009    0.316
## t           -1.662e+02  1.619e+02  -1.026    0.308
## t2           8.577e-02  8.215e-02   1.044    0.299
## t3          -1.476e-05  1.389e-05  -1.063    0.291
## 
## Residual standard error: 1.814 on 86 degrees of freedom
## Multiple R-squared:  0.7425, Adjusted R-squared:  0.7335 
## F-statistic: 82.67 on 3 and 86 DF,  p-value: < 2.2e-16

From the above two Quadratic model we see that the R². value seems to be improved and time components are significant for Quadratic model of order 2. However this is not the same when we increase the oder of the quadratic model to order of 3. Hence we would say Quadratic model os order 2 seems to be better compared to the remaining models tried.

Harmonic Models

Let us try modelling an harmonic model eventhough there is no apparent seasonality effect considering it is a yearly data points. We can still explore and see how the model is performing for the given data points:-

har.=harmonic(ozoneLayer,0.45) # calculate cos(2*pi*t) and sin(2*pi*t)
harmonicModel=lm(ozoneLayer~har.)

## 
## Call:
## lm(formula = ozoneLayer ~ har.)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3520 -1.8905  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.4441

#har.=harmonic(ozoneLayer,0.45) # calculate cos(2*pi*t) and sin(2*pi*t)
harmonicModel=lm(ozoneLayer~har.+t)

## 
## Call:
## lm(formula = ozoneLayer ~ har. + t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5565 -1.5892  0.0691  1.5110  4.7170 
## 
## Coefficients: (1 not defined because of singularities)
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.133e+02  1.633e+01  13.057   <2e-16 ***
## har.cos(2*pi*t)         NA         NA      NA       NA    
## har.sin(2*pi*t)  2.410e+11  4.122e+11   0.585     0.56    
## t               -1.098e-01  8.289e-03 -13.241   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.04 on 87 degrees of freedom
## Multiple R-squared:  0.6706, Adjusted R-squared:  0.663 
## F-statistic: 88.54 on 2 and 87 DF,  p-value: < 2.2e-16

The harmonic models tried gave a different results with the first harmonic model being insignificant with p value > 0.05 whereas the second harmonic model is significant considering the p-value. But the point to be noted here is that the second harmonic model does not have either the cosine or the sine components significant. It is better to go with wither the linear or the quadratic model as they both give us a good fit compared to the rest.

Model Validation

Residual Analysis

We can perform a Residual Analysis to check the model validation.

Time Series Plot

par(mar=c(2,4,2,1),mfrow=c(2,2),cex.main=0.75, cex.lab=0.75, cex.axis=0.75)
plot(y=rstudent(lmModel),x=as.vector(time(ozoneLayer)), xlab='Time',
     ylab='Standardized Residuals',type='o', 
     main = "Time series plot of residuals of Linear Model")
plot(y=rstudent(quadModel_2),x=as.vector(time(ozoneLayer)), xlab='Time',
     ylab='Standardized Residuals ',type='o', 
     main = "Time series plot of residuals of Quadratic Model of Order 2")
plot(y=rstudent(quadModel_3),x=as.vector(time(ozoneLayer)), xlab='Time',
     ylab='Standardized Residuals ',type='o', 
     main = "Time series plot of residuals of Quadratic Model of Order 3")
plot(y=rstudent(harmonicModel),x=as.vector(time(ozoneLayer)), xlab='Time',
     ylab='Standardized Residuals ',type='o', 
     main = "Time series plot of residuals of Harmonic Model of Order 2")

Time Series plot of the Standardised Residuals of the models explored

From Fig 5. we see that all the residuals are hovering over the mean levels(mean 0), for all the above series but we could look into the normal distribution of the residuals using the below Histograms.

Histograms of the standardised residuals

par(mar=c(2,4,2,1),mfrow=c(2,2),cex.main=0.75, cex.lab=0.75, cex.axis=0.75)
hist(rstudent(lmModel),xlab='Standardized Residuals for Linear Models', freq = FALSE, ylim = c(0,0.6)
     ,main='Standardized Residuals for Linear Models')
curve(dnorm(x, mean=mean(rstudent(lmModel)), sd=sd(rstudent(lmModel))), col="red", lwd=2, add=TRUE, yaxt="n")

hist(rstudent(quadModel_2),xlab='Standardized Residuals for Quadratic Models', freq = FALSE, ylim = c(0,0.6)
     ,main='Standardized Residuals for Quadratic Model of Order 2')
curve(dnorm(x, mean=mean(rstudent(quadModel_2)), sd=sd(rstudent(quadModel_2))), col="red", lwd=2, add=TRUE, yaxt="n")

hist(rstudent(quadModel_3),xlab='Standardized Residuals for Quadratic Models', freq = FALSE, ylim = c(0,0.6)
     ,main='Standardized Residuals for Quadratic Model of Order 3')
curve(dnorm(x, mean=mean(rstudent(quadModel_3)), sd=sd(rstudent(quadModel_3))), col="red", lwd=2, add=TRUE, yaxt="n")

hist(rstudent(harmonicModel),xlab='Standardized Residuals for Quadratic Models', freq = FALSE, ylim = c(0,0.6)
     ,main='Standardized Residuals for Harmonic Model')
curve(dnorm(x, mean=mean(rstudent(harmonicModel)), sd=sd(rstudent(harmonicModel))), col="red", lwd=2, add=TRUE, yaxt="n")

Histograms of the Standardised Residuals of the models explored

From above Fig 6., it is clear that the Linear and the Quadratic model of order 2 are apporximately normal compared to the remaining two, eventhough we notice soem outliers in these models.

Normal QQ plots of the standardised residuals

par(mar=c(2,4,2,1),mfrow=c(2,2),cex.main=0.75, cex.lab=0.75, cex.axis=0.75)

qqnorm(rstudent(lmModel), main="QQ plot for linear Model")
qqline(rstudent(lmModel), col = 2, lwd = 1, lty = 2)

qqnorm(rstudent(quadModel_2), main="QQ plot for Quadratic Model of Order 2")
qqline(rstudent(quadModel_2), col = 2, lwd = 1, lty = 2)

qqnorm(rstudent(quadModel_3), main="QQ plot for Quadratic Model of Order 3")
qqline(rstudent(quadModel_3), col = 2, lwd = 1, lty = 2)

qqnorm(rstudent(harmonicModel), main="QQ plot for Harmonic Model")
qqline(rstudent(harmonicModel), col = 2, lwd = 1, lty = 2)

Normal QQ plots of the Standardised Residuals of the models explored

From Fig 7., it is interesting to see that only the Quadratic model of order 2 seems to be having a approximately normal fits, followed by the linear model.

Shapiro Wilk - Test

shapiro.test(rstudent(lmModel))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(lmModel)
## W = 0.98733, p-value = 0.5372

shapiro.test(rstudent(quadModel_2))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(quadModel_2)
## W = 0.98889, p-value = 0.6493

shapiro.test(rstudent(quadModel_3))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(quadModel_3)
## W = 0.98966, p-value = 0.7064

shapiro.test(rstudent(harmonicModel))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(harmonicModel)
## W = 0.98546, p-value = 0.4166

The above tests prove evident that there is no enough evidence to reject the null Hypothesis, hence the residuals are normal.

ACF Plots

Finally, we can plot a ACF to validate the significant correlation in the residuals.

par(mar=c(2,4,2,1),mfrow=c(2,2),cex.main=0.75, cex.lab=0.75, cex.axis=0.75)
acf(lmModel$residuals, main = "ACF of standardized residuals of linear data")
acf(quadModel_2$residuals, main = "ACF of standardized residuals of Quadratic model of Order 2")
acf(quadModel_3$residuals, main = "ACF of standardized residuals of Quadratic model of Order 3")
acf(harmonicModel$residuals, main = "ACF of standardized residuals of Harmonic model")

ACF plots of the Standardised Residuals of the models explored

From the ACF plots Fig 8. we see that there is a significant correlation in certain lags of residuals for linear and quadratic models with the auto-correlation being slightly higher for the Qudratic plots.

Interpretation from our Analysis

From our analysis we believe that the Quadratic Model of order 2 has a good fit followed by the Linear Model. But considering the overall statistics results i.e. considering the residual analysis we cannot say that this is a robust model. We might have to explore more complex models and consider performng some transformations to have a robust best model. However the scope of this assignment is to perform a trend model we ;eave that as a limitation & future scope of work.

Forecast - 5 step ahead

In this final step we will predict a 5 year ahead forecast of the Ozone Layer Thickness changes with two of our good models Linear Model & Quadratic Model of Order 2.

h = 5 # 5 steps ahed forecasts
ta = seq(end(ozoneLayer)[1], (end(ozoneLayer)[1]+h), 1)
#new = data.frame(ta)
forecasts_lm = predict(lmModel, data.frame(ta), interval = "prediction")

ta2 = ta^2
#newQuad = data.frame(ta ,ta2)
forecasts_quad = predict(quadModel_2, data.frame(ta ,ta2), interval = "prediction")
#print(head(forecasts_lm))

Predicted changes in Ozone Layer Thickness in Dobston Units
fit	lwr	upr
1.6940317	-2.431767	5.8198300
1.5840026	-2.538931	5.7069365
1.4739734	-2.646159	5.5941061
1.3639443	-2.753450	5.4813390
1.2539152	-2.860805	5.3686354
1.1438861	-2.968223	5.2559953
1.0338569	-3.075705	5.1434189
0.9238278	-3.183251	5.0309063
0.8137987	-3.290860	4.9184576
0.7037696	-3.398534	4.8060729
0.5937404	-3.506271	4.6937523
0.4837113	-3.614073	4.5814960
0.3736822	-3.721940	4.4693041
0.2636531	-3.829871	4.3571766
0.1536239	-3.937866	4.2451136
0.0435948	-4.045926	4.1331153
-0.0664343	-4.154050	4.0211817
-0.1764634	-4.262240	3.9093130
-0.2864926	-4.370494	3.7975092
-0.3965217	-4.478814	3.6857703
-0.5065508	-4.587198	3.5740966
-0.6165799	-4.695648	3.4624880
-0.7266091	-4.804163	3.3509447
-0.8366382	-4.912743	3.2394666
-0.9466673	-5.021389	3.1280540
-1.0566964	-5.130100	3.0167067
-1.1667256	-5.238876	2.9054250
-1.2767547	-5.347718	2.7942089
-1.3867838	-5.456626	2.6830584
-1.4968129	-5.565599	2.5719736
-1.6068421	-5.674639	2.4609545
-1.7168712	-5.783743	2.3500012
-1.8269003	-5.892914	2.2391137
-1.9369294	-6.002151	2.1282920
-2.0469586	-6.111453	2.0175363
-2.1569877	-6.220822	1.9068466
-2.2670168	-6.330256	1.7962228
-2.3770459	-6.439757	1.6856650
-2.4870751	-6.549323	1.5751732
-2.5971042	-6.658956	1.4647475
-2.7071333	-6.768654	1.3543878
-2.8171624	-6.878419	1.2440942
-2.9271916	-6.988250	1.1338668
-3.0372207	-7.098147	1.0237054
-3.1472498	-7.208110	0.9136102
-3.2572789	-7.318139	0.8035811
-3.3673081	-7.428234	0.6936180
-3.4773372	-7.538396	0.5837212
-3.5873663	-7.648623	0.4738904
-3.6973954	-7.758917	0.3641257
-3.8074246	-7.869276	0.2544271
-3.9174537	-7.979702	0.1447946
-4.0274828	-8.090194	0.0352281
-4.1375119	-8.200752	-0.0742724
-4.2475411	-8.311375	-0.1837068
-4.3575702	-8.422065	-0.2930753
-4.4675993	-8.532821	-0.4023778
-4.5776284	-8.643642	-0.5116145
-4.6876576	-8.754530	-0.6207852
-4.7976867	-8.865483	-0.7298901
-4.9077158	-8.976502	-0.8389293
-5.0177449	-9.087587	-0.9479027
-5.1277741	-9.198738	-1.0568105
-5.2378032	-9.309954	-1.1656526
-5.3478323	-9.421235	-1.2744291
-5.4578614	-9.532583	-1.3831402
-5.5678906	-9.643995	-1.4917857
-5.6779197	-9.755473	-1.6003660
-5.7879488	-9.867017	-1.7088809
-5.8979779	-9.978625	-1.8173305
-6.0080071	-10.090299	-1.9257150
-6.1180362	-10.202038	-2.0340344
-6.2280653	-10.313842	-2.1422889
-6.3380944	-10.425710	-2.2504784
-6.4481235	-10.537644	-2.3586030
-6.5581527	-10.649642	-2.4666630
-6.6681818	-10.761705	-2.5746583
-6.7782109	-10.873833	-2.6825890
-6.8882400	-10.986025	-2.7904553
-6.9982692	-11.098281	-2.8982573
-7.1082983	-11.210602	-3.0059950
-7.2183274	-11.322986	-3.1136685
-7.3283565	-11.435435	-3.2212781
-7.4383857	-11.547948	-3.3288237
-7.5484148	-11.660524	-3.4363056
-7.6584439	-11.773164	-3.5437237
-7.7684730	-11.885868	-3.6510783
-7.8785022	-11.998635	-3.7583695
-7.9885313	-12.111465	-3.8655974
-8.0985604	-12.224359	-3.9727621

Predicted changes in Ozone Layer Thickness in Dobston Units
fit	lwr	upr
-0.3035205	-4.079911	3.4728701
-0.2788832	-4.041082	3.4833152
-0.2573065	-4.006376	3.4917627
-0.2387904	-3.975750	3.4981690
-0.2233349	-3.949160	3.5024903
-0.2109400	-3.926563	3.5046829
-0.2016057	-3.907914	3.5047030
-0.1953319	-3.893171	3.5025070
-0.1921188	-3.882290	3.4980519
-0.1919663	-3.875227	3.4912946
-0.1948744	-3.871942	3.4821929
-0.2008431	-3.872391	3.4707048
-0.2098724	-3.876534	3.4567894
-0.2219623	-3.884331	3.4404062
-0.2371128	-3.895741	3.4215157
-0.2553239	-3.910727	3.4000790
-0.2765956	-3.929250	3.3760586
-0.3009279	-3.951273	3.3494175
-0.3283208	-3.976762	3.3201201
-0.3587743	-4.005680	3.2881317
-0.3922884	-4.037996	3.2534189
-0.4288631	-4.073675	3.2159492
-0.4684984	-4.112688	3.1756915
-0.5111943	-4.155005	3.1326160
-0.5569508	-4.200595	3.0866938
-0.6057679	-4.249433	3.0378976
-0.6576456	-4.301493	2.9862013
-0.7125839	-4.356748	2.9315799
-0.7705828	-4.415176	2.8740098
-0.8316423	-4.476754	2.8134690
-0.8957624	-4.541461	2.7499363
-0.9629432	-4.609279	2.6833922
-1.0331845	-4.680187	2.6138183
-1.1064864	-4.754170	2.5411977
-1.1828489	-4.831213	2.4655147
-1.2622720	-4.911299	2.3867549
-1.3447557	-4.994417	2.3049052
-1.4303000	-5.080554	2.2199540
-1.5189050	-5.169701	2.1318908
-1.6105705	-5.261848	2.0407065
-1.7052966	-5.356986	1.9463933
-1.8030833	-5.455111	1.8489447
-1.9039306	-5.556217	1.7483556
-2.0078386	-5.660299	1.6446219
-2.1148071	-5.767355	1.5377412
-2.2248362	-5.877384	1.4277120
-2.3379259	-5.990386	1.3145345
-2.4540763	-6.106363	1.1982099
-2.5732872	-6.225315	1.0787409
-2.6955587	-6.347249	0.9561312
-2.8208909	-6.472168	0.8303861
-2.9492836	-6.600079	0.7015122
-3.0807369	-6.730991	0.5695171
-3.2152508	-6.864912	0.4344101
-3.3528254	-7.001852	0.2962015
-3.4934605	-7.141824	0.1549031
-3.6371562	-7.284840	0.0105279
-3.7839126	-7.430915	-0.1369098
-3.9337295	-7.580065	-0.2873942
-4.0866071	-7.732306	-0.4409083
-4.2425452	-7.887657	-0.5974339
-4.4015439	-8.046137	-0.7569513
-4.5636033	-8.207767	-0.9194395
-4.7287232	-8.372570	-1.0848764
-4.8969038	-8.540569	-1.2532382
-5.0681449	-8.711790	-1.4245003
-5.2424467	-8.886257	-1.5986364
-5.4198090	-9.063999	-1.7756191
-5.6002319	-9.245044	-1.9554197
-5.7837155	-9.429423	-2.1380082
-5.9702596	-9.617166	-2.3233536
-6.1598644	-9.808305	-2.5114235
-6.3525297	-10.002875	-2.7021843
-6.5482557	-10.200910	-2.8956015
-6.7470422	-10.402445	-3.0916393
-6.9488894	-10.607518	-3.2902609
-7.1537971	-10.816166	-3.4914286
-7.3617655	-11.028427	-3.6951037
-7.5727945	-11.244342	-3.9012465
-7.7868840	-11.463951	-4.1098168
-8.0040342	-11.687295	-4.3207733
-8.2242449	-11.914416	-4.5340742
-8.4475163	-12.145355	-4.7496773
-8.6738483	-12.380157	-4.9675396
-8.9032408	-12.618864	-5.1876180
-9.1356940	-12.861519	-5.4098688
-9.3712077	-13.108167	-5.6342484
-9.6097821	-13.358851	-5.8607129
-9.8514171	-13.613615	-6.0892187
-10.0961126	-13.872503	-6.3197220

plot(ozoneLayer, ylab = "Ozone Layer Thickness in Dobson Units)",ylim=c(-15,5),type="o",
     main = " Time Series Plot of Ozone Layer Thickness - 5 Year Forecasts")

lines(ts(as.vector(forecasts_lm[,1]), start = 2017), col="red", type="l", lwd=3)
lines(ts(as.vector(forecasts_lm[,2]), start = 2017), col="blue", type="l", lwd=3)
lines(ts(as.vector(forecasts_lm[,3]), start = 2017), col="blue", type="l", lwd=3)

legend("bottomleft", lty=1, pch=1, 
       col=c("black","red", "blue"), text.width = 18,
       c("OzoneData","LinearModel - Fitted forecast", "LinearModel - 5% forecast limits"), 
       bty = "n")

Plotting of Time Series data and the 5 year ahead forecasts[Linear Model]

plot(ozoneLayer, ylab = "Ozone Layer Thickness in Dobson Units)",ylim=c(-15,5),type="o",
     main = " Time Series Plot of Ozone Layer Thickness - 5 Year Forecasts")

lines(ts(as.vector(forecasts_quad[,1]), start = 2017), col="magenta", type="l", lwd=3)
lines(ts(as.vector(forecasts_quad[,2]), start = 2017), col="brown", type="l", lwd=3)
lines(ts(as.vector(forecasts_quad[,3]), start = 2017), col="brown", type="l", lwd=3)


legend("bottomleft", lty=1, pch=1, 
       col=c("black","magenta","brown"), text.width = 18,
       c("OzoneData", "QuadraticModel - Fitted forecast", "QuadraticModel - 5% forecast limits"), 
       bty = "n")

Plotting of Time Series data and the 5 year ahead forecasts[Quadratic Model]

Conclusion

From the above Fig 9. & Fig 10. we observe that the Quadratic Model gives a close enough prediction visually compared to the Linear Model with prediction looking like an intervention point which may not be possible. However, we cannot rule out the option and hence mentioning it as a future scope of work.

References

Rob J Hyndman. 2018. Time Series. [ONLINE] Available at: https://robjhyndman.com/categories/time-series/. [Accessed 31 March 2018].
H., R., 2017. Time Series Analysis and Its Applications: With R Examples (Springer Texts in Statistics). Springer.
D., J., 2009. Time Series Analysis. Springer-Verlag New York, LLC.

Ozone Layer Thickness Analysis

MATH1318 Time Series Analysis- Assignment 1

Arjun Balaji (s3629999)

31 March 2019

Problem

Data Description

Descriptive Analysis

Model

Linear Regression Model

Quadratic Models

Harmonic Models

Model Validation

Residual Analysis

Time Series Plot

Histograms of the standardised residuals

Normal QQ plots of the standardised residuals

Shapiro Wilk - Test

ACF Plots

Interpretation from our Analysis

Forecast - 5 step ahead

Conclusion

References