Final Project

1. Introduction

For this project we have chosen to analyse and forecast on the air quality of the city of Madrid, Spain. In the dataset there are many types of pollutants present in the air but we have chosen to focus on the levels of the ozone gas in the air. High ozone levels can lead to breating difficulty, lung damage and aggravate existing lung diseases. As such we feel this is an important variable to monitor. To determine the severity of the future damage caused by this gas we will be forecasting the levels for the next 10 months.

2. The Data

The dataset used in this forecasting report is ’Air Quality of Madrid (2001 - 2018), it contains several measurements for air pollutants, among them the one variable we are interested in, the ozone level.

Source: https://www.kaggle.com/decide-soluciones/air-quality-madrid

3. Packages

Below are the packages that will be used to conduct this analysis.

library(readr)
library(x12)
library(dynlm)
library(dLagM)
library(forecast)
library(AER)
library(ggplot2)
library(TSA)
library(Hmisc)
library(tseries)
library(dplyr)

4. Preparing the Data

First we need to load up the 18 datasets that will be used. Each dataset corresponds to a year between 2001 and 2018.

year2001 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2001.csv")
year2002 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2002.csv")
year2003 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2003.csv")
year2004 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2004.csv")
year2005 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2005.csv")
year2006 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2006.csv")
year2007 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2007.csv")
year2008 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2008.csv")
year2009 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2009.csv")
year2010 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2010.csv")
year2011 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2011.csv")
year2012 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2012.csv")
year2013 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2013.csv")
year2014 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2014.csv")
year2015 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2015.csv")
year2016 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2016.csv")
year2017 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2017.csv")
year2018 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2018.csv")
stations <- read_csv("air-quality-madrid/stations.csv")

Next we need to modify the format of the date column in order to use them as time series data later.

Then we combine all 18 datasets using the rbind() function to generate the final dataset we will use for our analysis and forecasting. Preview of the final dataset below:

5. Visualisation

In this section we will begin by visualising the daily ozone levels for Madrid in Figure 1 to get a general idea of how the levels vary.

For the sake of predicting the ozone levels 10 months ahead, we will be converting the values to a monthly figure, the mean ozone levels for the month, and visualising again to better suit the goal of this project. Figure 2 provides the monthly ozone levels at 6 month intervals for the years 2001 to 2018.

6. Time Series Analysis

With basic visualisations complete, time series analysis can now be conducted on the final dataset. First, the dataset is converted to a time series and then an initial visualisation (Figure 3) is done to look for any obvious characteristics.

From figure 3 an obvious seasonal trend can be observed based on the behaviour of the points throughout the years. It can also be observed that the months July and August contain the highest ozone levels whereas December has the lowest levels over the time period. No notable intervention points observed. There is a changing variance with the mean, along with a large fluctuation and auto-regressive hehaviour.

6.1 Checking for Stationarity

Using an ADF test we shall determine if the dataset is stationary or not.

As the above result shows, the p value indicates this time series is stationary.

6.2 ACF, PACF & First Decomposition

From figure 3, there is already an observed seasonality for the dataset but just to be sure, an ACF and PACF plot will be made to check (Figure 4 & 5).

From the above figures we can see that there is in fact seasonality within the dataset.

Next we apply the x12 Decomposition to observe the decomposed lines.

There is auto-regressive behaviour which fluctuates around the mean and is present in the figure acroos the entire time series.

From figure 6, the seasonally adjusted keeps along with the trend and orginal trend but with some flucuations, which indicates the seasonality is not completey removed. Therefore, extra differencing is required. Manwhile, it is hard to conclude the exsistence of multi-seasonality beaviour. Furthermore, there is evidence the trend slighlty increases with time. The trend effect will be investigated in the next part.

6.3 Further Decomposition

First an STL decomposition was run on the time series to observe the factors of the time series seperately (Figure 7). Then in Figure 8, the seasonal component is extracted and plotted, showing a much more obvious trend pattern emerging in the time series. Following Figure 8, the trend component is removed as well giving us Figure 9 in which we can see that the fluctuations around the mean are more equal throughout the years and using an ADF test, the time series is concluded to be stationary.

7. Model Fitting

Model fitting will be divided into multiple sections with each section corresponding to specific types of models. Results from the MASE and AIC will be summarised in a table at the bottom in order to determine which model is the most accurate amongst the different types of models.

7.1 Dynamic Linear Models

Model Name	AIC	MASE
\(model.dynlm1\)	1231.498	0.4039228
\(model.dynlm1.2\)	1263.965	0.4313237
\(model.dynlm1.3\)	1561.645	0.9858214
\(model.dynlm2\)	1255.370	0.4322566
\(model.dynlm4\)	1245.226	0.4173508

From the table and based on the AIC, BIC and MASE results, model.dynlm1 is the best model as it has the lowest AIC, BIC and MASE.

However, from the residuals, there is one significant value in the acf and the tails of the residuals histogram are a bit longer than we would prefer and it shows to be random in the line chart in the residuals check.

7.2 Exponential Smoothing Models

Model Name	AIC	MASE
\(model.simple\)	2112.865	1.649948
\(model.hw1\)	1778.609	0.687281
\(model.hw2\)	1779.762	0.684349
\(model.hw3\)	1814.686	0.697647
\(model.hw4\)	1829.954	0.707525

From the table above, we can see that the best model is model.hw2 which is the Damped Holt-Winters’ additive method with a MASE value of 0.6843495. Looking at the residuals as well, the histogram appears mostly normal with short tails, there are only 2 significant values in the ACF and the line chart shows equal variances around the mean.

7.3 State-Space Models

Model Name	AIC	MASE
\(model.ANA\)	1775.858	0.6847280
\(model.AAA\)	1778.488	0.6872814
\(model.AAdA.damped\)	1779.938	0.6843495
\(model.ANN\)	2112.889	1.6495840
\(model.MAM\)	1815.027	0.6918962
\(model.MAM.damped\)	1815.027	0.6918962
\(model.MMM\)	1814.531	0.6919849
\(model.MMdM\)	1814.531	0.6919849
\(model.ZZZ\)	1775.858	0.6847280

From the table above the best model is model.ANA which uses the state space ANA model and has the lowest AIC, BIC and MASE values among the other models. Looking at the residuals, there are 2 significant values in the ACF, the histograms appears mostly normal but with some long tails and the line plot has an almost constant variance around the mean.

8. Model selection for Forecasting

Taking the best models from each modelling category, we will now compare them in the table below.

Model Name	AIC	MASE
\(model.dynlm1\)	1231.498	0.4039228
\(model.hw2\)	1779.762	0.684349
\(model.ANA\)	1775.858	0.6847280

Finally, from the table above comparing the best models from each category, we have now determined that the best model to use for forecasting is the dynamic linear model, model.dynlm1.

9. Forecasting

Below is the 10 Month ahead forecast of Ozone levels in the City of Madrid using the selected dynamic linear model.

10. Conclusion

To conclude, out of the three modelling techniques applied, our best model resulted from a dynamic linear model with a lag of 1 and the trend and seasonal components added as it resulted in the lowest AIC and MASE values out of all the models.

Based on the forecast above, we are expecting to see an increase in ozone levels for the city of Madrid in the next 10 months, as expected, and residents should take steps to guard themselves against the potential hazards the gas presents and potentially make efforts to enact changes to reduce ozone gas levels for the future.

11. Code

#3
library(readr)
library(x12)
library(dynlm)
library(dLagM)
library(forecast)
library(AER)
library(ggplot2)
library(TSA)
library(Hmisc)
library(tseries)
library(dplyr)

#4
year2001 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2001.csv")
year2002 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2002.csv")
year2003 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2003.csv")
year2004 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2004.csv")
year2005 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2005.csv")
year2006 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2006.csv")
year2007 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2007.csv")
year2008 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2008.csv")
year2009 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2009.csv")
year2010 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2010.csv")
year2011 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2011.csv")
year2012 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2012.csv")
year2013 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2013.csv")
year2014 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2014.csv")
year2015 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2015.csv")
year2016 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2016.csv")
year2017 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2017.csv")
year2018 <- read_csv("air-quality-madrid/csvs_per_year/madrid_2018.csv")
stations <- read_csv("air-quality-madrid/stations.csv")

##Change date format and selecting the O3 variables (pollutants)
##Since the date format is in **factor type**, I change those type into Posixct and selecting only O3 variables.

year2001$date <- as.POSIXct(year2001$date, format = "%Y-%m-%d %H:%M:%S")
year2001$Date <- format(year2001$date, "%Y-%m-%d")
year2001$Time <- format(year2001$date, "%T")
year2001$Date <- as.POSIXct(year2001$Date, format = "%Y-%m-%d")

year2001$day <- format(year2001$date, "%d")
year2001$month <- format(year2001$date, "%m")
year2001$year <- format(year2001$date, "%Y")
year2001$hour <- format(year2001$date, "%H")

m2001 <- year2001 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2002$date <- as.POSIXct(year2002$date, format = "%Y-%m-%d %H:%M:%S")
year2002$Date <- format(year2002$date, "%Y-%m-%d")
year2002$Time <- format(year2002$date, "%T")
year2002$Date <- as.POSIXct(year2002$Date, format = "%Y-%m-%d")

year2002$day <- format(year2002$date, "%d")
year2002$month <- format(year2002$date, "%m")
year2002$year <- format(year2002$date, "%Y")
year2002$hour <- format(year2002$date, "%H")

m2002 <- year2002 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2003$date <- as.POSIXct(year2003$date, format = "%Y-%m-%d %H:%M:%S")
year2003$Date <- format(year2003$date, "%Y-%m-%d")
year2003$Time <- format(year2003$date, "%T")
year2003$Date <- as.POSIXct(year2003$Date, format = "%Y-%m-%d")

year2003$day <- format(year2003$date, "%d")
year2003$month <- format(year2003$date, "%m")
year2003$year <- format(year2003$date, "%Y")
year2003$hour <- format(year2003$date, "%H")

m2003 <- year2003 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2004$date <- as.POSIXct(year2004$date, format = "%Y-%m-%d %H:%M:%S")
year2004$Date <- format(year2004$date, "%Y-%m-%d")
year2004$Time <- format(year2004$date, "%T")
year2004$Date <- as.POSIXct(year2004$Date, format = "%Y-%m-%d")

year2004$day <- format(year2004$date, "%d")
year2004$month <- format(year2004$date, "%m")
year2004$year <- format(year2004$date, "%Y")
year2004$hour <- format(year2004$date, "%H")

m2004 <- year2004 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2005$date <- as.POSIXct(year2005$date, format = "%Y-%m-%d %H:%M:%S")
year2005$Date <- format(year2005$date, "%Y-%m-%d")
year2005$Time <- format(year2005$date, "%T")
year2005$Date <- as.POSIXct(year2005$Date, format = "%Y-%m-%d")

year2005$day <- format(year2005$date, "%d")
year2005$month <- format(year2005$date, "%m")
year2005$year <- format(year2005$date, "%Y")
year2005$hour <- format(year2005$date, "%H")

m2005 <- year2005 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2006$date <- as.POSIXct(year2006$date, format = "%Y-%m-%d %H:%M:%S")
year2006$Date <- format(year2006$date, "%Y-%m-%d")
year2006$Time <- format(year2006$date, "%T")
year2006$Date <- as.POSIXct(year2006$Date, format = "%Y-%m-%d")

year2006$day <- format(year2006$date, "%d")
year2006$month <- format(year2006$date, "%m")
year2006$year <- format(year2006$date, "%Y")
year2006$hour <- format(year2006$date, "%H")

m2006 <- year2006 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2007$date <- as.POSIXct(year2007$date, format = "%Y-%m-%d %H:%M:%S")
year2007$Date <- format(year2007$date, "%Y-%m-%d")
year2007$Time <- format(year2007$date, "%T")
year2007$Date <- as.POSIXct(year2007$Date, format = "%Y-%m-%d")

year2007$day <- format(year2007$date, "%d")
year2007$month <- format(year2007$date, "%m")
year2007$year <- format(year2007$date, "%Y")
year2007$hour <- format(year2007$date, "%H")

m2007 <- year2007 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2008$date <- as.POSIXct(year2008$date, format = "%Y-%m-%d %H:%M:%S")
year2008$Date <- format(year2008$date, "%Y-%m-%d")
year2008$Time <- format(year2008$date, "%T")
year2008$Date <- as.POSIXct(year2008$Date, format = "%Y-%m-%d")

year2008$day <- format(year2008$date, "%d")
year2008$month <- format(year2008$date, "%m")
year2008$year <- format(year2008$date, "%Y")
year2008$hour <- format(year2008$date, "%H")

m2008 <- year2008 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2009$date <- as.POSIXct(year2009$date, format = "%Y-%m-%d %H:%M:%S")
year2009$Date <- format(year2009$date, "%Y-%m-%d")
year2009$Time <- format(year2009$date, "%T")
year2009$Date <- as.POSIXct(year2009$Date, format = "%Y-%m-%d")

year2009$day <- format(year2009$date, "%d")
year2009$month <- format(year2009$date, "%m")
year2009$year <- format(year2009$date, "%Y")
year2009$hour <- format(year2009$date, "%H")

m2009 <- year2009 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2010$date <- as.POSIXct(year2010$date, format = "%Y-%m-%d %H:%M:%S")
year2010$Date <- format(year2010$date, "%Y-%m-%d")
year2010$Time <- format(year2010$date, "%T")
year2010$Date <- as.POSIXct(year2010$Date, format = "%Y-%m-%d")

year2010$day <- format(year2010$date, "%d")
year2010$month <- format(year2010$date, "%m")
year2010$year <- format(year2010$date, "%Y")
year2010$hour <- format(year2010$date, "%H")

m2010 <- year2010 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2011$date <- as.POSIXct(year2011$date, format = "%Y-%m-%d %H:%M:%S")
year2011$Date <- format(year2011$date, "%Y-%m-%d")
year2011$Time <- format(year2011$date, "%T")
year2011$Date <- as.POSIXct(year2011$Date, format = "%Y-%m-%d")

year2011$day <- format(year2011$date, "%d")
year2011$month <- format(year2011$date, "%m")
year2011$year <- format(year2011$date, "%Y")
year2011$hour <- format(year2011$date, "%H")

m2011 <- year2011 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2012$date <- as.POSIXct(year2012$date, format = "%Y-%m-%d %H:%M:%S")
year2012$Date <- format(year2012$date, "%Y-%m-%d")
year2012$Time <- format(year2012$date, "%T")
year2012$Date <- as.POSIXct(year2012$Date, format = "%Y-%m-%d")

year2012$day <- format(year2012$date, "%d")
year2012$month <- format(year2012$date, "%m")
year2012$year <- format(year2012$date, "%Y")
year2012$hour <- format(year2012$date, "%H")

m2012 <- year2012 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2013$date <- as.POSIXct(year2013$date, format = "%Y-%m-%d %H:%M:%S")
year2013$Date <- format(year2013$date, "%Y-%m-%d")
year2013$Time <- format(year2013$date, "%T")
year2013$Date <- as.POSIXct(year2013$Date, format = "%Y-%m-%d")

year2013$day <- format(year2013$date, "%d")
year2013$month <- format(year2013$date, "%m")
year2013$year <- format(year2013$date, "%Y")
year2013$hour <- format(year2013$date, "%H")

m2013 <- year2013 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2014$date <- as.POSIXct(year2014$date, format = "%Y-%m-%d %H:%M:%S")
year2014$Date <- format(year2014$date, "%Y-%m-%d")
year2014$Time <- format(year2014$date, "%T")
year2014$Date <- as.POSIXct(year2014$Date, format = "%Y-%m-%d")

year2014$day <- format(year2014$date, "%d")
year2014$month <- format(year2014$date, "%m")
year2014$year <- format(year2014$date, "%Y")
year2014$hour <- format(year2014$date, "%H")

m2014 <- year2014 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2015$date <- as.POSIXct(year2015$date, format = "%Y-%m-%d %H:%M:%S")
year2015$Date <- format(year2015$date, "%Y-%m-%d")
year2015$Time <- format(year2015$date, "%T")
year2015$Date <- as.POSIXct(year2015$Date, format = "%Y-%m-%d")

year2015$day <- format(year2015$date, "%d")
year2015$month <- format(year2015$date, "%m")
year2015$year <- format(year2015$date, "%Y")
year2015$hour <- format(year2015$date, "%H")

m2015 <- year2015 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2016$date <- as.POSIXct(year2016$date, format = "%Y-%m-%d %H:%M:%S")
year2016$Date <- format(year2016$date, "%Y-%m-%d")
year2016$Time <- format(year2016$date, "%T")
year2016$Date <- as.POSIXct(year2016$Date, format = "%Y-%m-%d")

year2016$day <- format(year2016$date, "%d")
year2016$month <- format(year2016$date, "%m")
year2016$year <- format(year2016$date, "%Y")
year2016$hour <- format(year2016$date, "%H")

m2016 <- year2016 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2017$date <- as.POSIXct(year2017$date, format = "%Y-%m-%d %H:%M:%S")
year2017$Date <- format(year2017$date, "%Y-%m-%d")
year2017$Time <- format(year2017$date, "%T")
year2017$Date <- as.POSIXct(year2017$Date, format = "%Y-%m-%d")

year2017$day <- format(year2017$date, "%d")
year2017$month <- format(year2017$date, "%m")
year2017$year <- format(year2017$date, "%Y")
year2017$hour <- format(year2017$date, "%H")

m2017 <- year2017 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)


year2018$date <- as.POSIXct(year2018$date, format = "%Y-%m-%d %H:%M:%S")
year2018$Date <- format(year2018$date, "%Y-%m-%d")
year2018$Time <- format(year2018$date, "%T")
year2018$Date <- as.POSIXct(year2018$Date, format = "%Y-%m-%d")

year2018$day <- format(year2018$date, "%d")
year2018$month <- format(year2018$date, "%m")
year2018$year <- format(year2018$date, "%Y")
year2018$hour <- format(year2018$date, "%H")

m2018 <- year2018 %>%
  select(date, Date, Time, day, month, year, hour, station, O_3)

madrid <- m2001 %>%
  rbind(m2002) %>%
  rbind(m2003) %>%
  rbind(m2004) %>%
  rbind(m2005) %>%
  rbind(m2006) %>%
  rbind(m2007) %>%
  rbind(m2008) %>%
  rbind(m2009) %>%
  rbind(m2010) %>%
  rbind(m2011) %>%
  rbind(m2012) %>%
  rbind(m2013) %>%
  rbind(m2014) %>%
  rbind(m2015) %>%
  rbind(m2016) %>%
  rbind(m2017) %>%
  rbind(m2018) 

radarstations <- stations %>% select(id, name)
madrid <-  right_join(madrid, radarstations, by = c("station" = "id"))
print(head(madrid))

#5
madrid_mean_daily <- madrid %>% group_by(Date) %>% summarise(O_3_mean = mean(O_3, na.rm = TRUE))
ggplot(madrid_mean_daily, aes(x = as.Date(Date), y = O_3_mean)) + geom_line(color = 'blue') +
  theme_bw() +
  labs(x = 'Year', y = 'Ozone Level (μg/m3)', title='Daily Ozone Levels in Madrid (Figure 1)') +
  theme(axis.text.x=element_text(size=10)) +
  scale_x_date(breaks = seq(as.Date("2001-01-01"), as.Date("2018-07-01"), by="12 months"), date_labels = "%Y")

madrid_mean_monthly <- madrid %>% group_by(year, month) %>% summarise( O_3_mean = mean(O_3, na.rm = TRUE))%>% mutate(time = paste(year, "-", month, "- 01"))

madrid_mean_monthly$time <- as.Date(madrid_mean_monthly$time, format = "%Y - %m - %d")

ggplot(madrid_mean_monthly, aes(x = time, y = O_3_mean)) + geom_line(color = 'blue') + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 90)) +
  labs(x = 'Month', y = 'Ozone Level (μg/m3)', title='Monthly Ozone Levels in Madrid (Figure 2)') +
  theme(axis.text.x=element_text(size=10)) +
  scale_x_date(breaks = seq(as.Date("2001-01-01"), as.Date("2018-07-01"), by="6 months"), date_labels = "%b-%y") 

#6
dummy1 <- as.vector(madrid_mean_monthly$O_3_mean)
madrid.ts <- ts(dummy1,start = c(2001,1), end = c(2018,5), frequency = 12)

plot(madrid.ts, ylab="Ozone Level (μg/m3)", main = "Time series plot for O3 (Figure 3)")
points(y=madrid.ts, x=time(madrid.ts), pch = as.vector(season(madrid.ts)))

#6.1
adf.test(madrid.ts)

#6.2
par(mfrow= c(1,2))
acf(madrid.ts, main = "Figure 4")
pacf(madrid.ts, main = "Figure 5")

#6.3
#Using seasonal differencing
fit.madrid <- stl(madrid.ts, t.window=15, s.window="periodic", robust=TRUE)
plot(fit.madrid, main = "Figure 7")
fit.madrid.seasonal = fit.madrid$time.series[,"seasonal"] 

#Extract the seansonal component from the output
madrid.seasonal.adjusted = madrid.ts - fit.madrid.seasonal
plot(madrid.seasonal.adjusted,xlab='Time', ylab ='Ozone Levels (μg/m3)', main = "Seasonal adjusted time series (Figure 8)")
points(y=madrid.seasonal.adjusted,x=time(madrid.seasonal.adjusted), pch=as.vector(season(madrid.seasonal.adjusted)))

#Extract the trend component from the output
fit.madrid.trend = fit.madrid$time.series[,"trend"] 
madrid.seasonal.trend.adjusted = madrid.ts - fit.madrid.seasonal - fit.madrid.trend
plot(madrid.seasonal.trend.adjusted,xlab='Time', ylab ='Ozone Level (μg/m3)', main = "Trend & seasonal adjusted time series (Figure 9)")
points(y=madrid.seasonal.trend.adjusted,x=time(madrid.seasonal.trend.adjusted), pch=as.vector(season(madrid.seasonal.trend.adjusted)))

adf.test(madrid.seasonal.trend.adjusted)

#7.1
#Dynlm Models
Y.t = madrid.ts

model.dynlm1 = dynlm(Y.t ~ L(Y.t , k = 1 ) + trend(Y.t) + season(Y.t))

model.dynlm2 = dynlm(Y.t ~ L(Y.t , k = 2 ) + trend(Y.t) + season(Y.t))

model.dynlm3 = dynlm(Y.t ~ L(Y.t , k = 1 ) + trend(Y.t) + season(Y.t)) #MASE 0.403897   

model.dynlm1.2 = dynlm(Y.t ~ L(Y.t , k = 1 ) + season(Y.t)) 

model.dynlm1.3 = dynlm(Y.t ~ L(Y.t , k = 1 ) + trend(Y.t))

model.dynlm4.1 = dynlm(Y.t ~ L(Y.t , k = 1 ) + L(Y.t , k = 2 )+ season(Y.t))

model.dynlm4.2 = dynlm(Y.t ~ L(Y.t , k = 1 ) + L(Y.t , k = 2 ) + season(Y.t))

#Creating table summary of model results
aic = AIC(model.dynlm1, model.dynlm1.2, model.dynlm1.3,model.dynlm2,model.dynlm3,model.dynlm4.1,model.dynlm4.2)
bic = BIC(model.dynlm1, model.dynlm1.2, model.dynlm1.3,model.dynlm2,model.dynlm3,model.dynlm4.1,model.dynlm4.2)
mase = MASE(lm(model.dynlm1), lm(model.dynlm1.2), lm(model.dynlm1.3),lm(model.dynlm2),lm(model.dynlm3),lm(model.dynlm4.1),lm(model.dynlm4.2))

dynlm_table <- cbind(aic,bic,mase)
dynlm_table
#Best model is model.dynlm1
summary(model.dynlm1)
checkresiduals(model.dynlm1)

#7.2
model.simple <- ses(madrid.ts,h=2*frequency(madrid.ts))#MASE 1.649948

model.hw1 <- hw(madrid.ts,seasonal="additive") #MASE 0.6872814 

model.hw2 <- hw(madrid.ts,seasonal="additive",damped = TRUE) #MASE 0.6843495 

model.hw3 <- hw(madrid.ts,seasonal="multiplicative") #MASE 0.6976475 

model.hw4 <- hw(madrid.ts,seasonal="multiplicative",exponential = TRUE) #MASE 0.7075254 

#Table creation
name_es <- c("model.simple","model.hw1","model.hw2","model.hw3","model.hw4")
mase_es <- c(1.649948,0.687281,0.684349,0.697647,0.707525)

table_es <- as.data.frame(cbind(name_es,mase_es))
table_es
#Best model is model.hw2
summary(model.hw2)
checkresiduals(model.hw2)

#7.3
model.ANA = ets(madrid.ts,model = "ANA") #MASE 0.684728

model.AAA = ets(madrid.ts,model = "AAA")#MASE 0.6872814

model.AAdA.damped = ets(madrid.ts,model = "AAA",damped = TRUE) #MASE 0.6843495

model.ANN = ets(madrid.ts,model = "ANN")#MASE 1.649584 

model.MAM = ets(madrid.ts,model = "MAM")#MASE 0.6918962

model.MAM.damped = ets(madrid.ts,model = "MAM",damped = TRUE) #MASE 0.6918962

model.MMM = ets(madrid.ts,model = "MMM")#MASE 0.6919849

model.MMdM = ets(madrid.ts,model = "MMM",damped = TRUE) #MASE 0.6919849

model.ZZZ = ets(madrid.ts, model="ZZZ")##MASE 0.684728

#Table creation
aic_ss <- AIC(model.ANA, model.AAA, model.AAdA.damped, model.ANN, model.MAM, model.MAM.damped, model.MMM, model.MMdM, model.ZZZ)
bic_ss <- BIC(model.ANA, model.AAA, model.AAdA.damped, model.ANN, model.MAM, model.MAM.damped, model.MMM, model.MMdM, model.ZZZ)
mase_ss <- c(0.684728,0.6872814,0.6843495,1.649584,0.6918962,0.6918962,0.6919849,0.6919849,0.684728)
table_ss <- as.data.frame(cbind(aic_ss,bic_ss,mase_ss))
table_ss
#Best model is model.ANA
summary(model.ANA)
checkresiduals(model.ANA)

#8
#Table creation for final model selection

model_final <- c("model.dynlm1", "model.hw2", "model.ANA")
type_final <- c("Dynamic Linear", "Exponential Smoothing", "State-Space")
aic_final <- c(1231.498, "NA", 1775.858)
bic_final <- c(1281.561, "NA", 1825.993)
mase_final <- c(0.4039228, 0.6843495, 0.6847280)

table_final <- as.data.frame(cbind(model_final,type_final,aic_final,bic_final,mase_final))
table_final

#9
#Dynlm Forecast
q = 10
n = nrow(model.dynlm1$model)

madrid.frc <- array(NA, (n+q))
madrid.frc[1:n] <- Y.t[1:length(Y.t)]

trend <- array(NA,q)
trend.start <- model.dynlm1$model[n, "trend(Y.t)"]
trend <- seq(trend.start, trend.start + q/12, 1/12)

Pi <- array(NA, dim = c(n+q,2))

for (i in 1:q){
  months <- array(0,11)
  months[(i-1)%%12] = 1
  data.new = c(1,madrid.frc[n-1+i],trend[i],months)
  madrid.frc[n+i] = as.vector(model.dynlm1$coefficients) %*% data.new
}

par(mfrow = c(1,1))
plot(Y.t,xlim=c(2001,2020), ylim=c(0,100),ylab='Ozone Level (μg/m3)',xlab='Year',main = "Dynlm Forecast of Ozone Levels in Madrid (Figure 10)")
lines(ts(madrid.frc[(n+1):(n+q)],start=c(2018,6),frequency = 12),col="green")
lines(ts(madrid.frc[(n+1):(n+q)] + 2*sd(madrid.frc[(n+1):(n+q)]), start = c(2018,6), frequency = 12), col = "red")
lines(ts(madrid.frc[(n+1):(n+q)] - 2*sd(madrid.frc[(n+1):(n+q)]), start = c(2018,6), frequency = 12), col = "blue")
legend("topleft", lty=1, pch=1, col=c("black","blue","red","green"), text.width = 4,
       c("Data","5% lower limit","95% upper limit","Mean prediction"))