RSProject2new.Rmd

<<<<<<< HEAD
---
title: "Project 2"
author: "Daniel, Tom, Rachael"
date: "Date"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
require("knitr")
sourcedir <- "C:/Users/Rachael/OneDrive - University of Virginia/UVA/Fall_2020/SYS_6021/Source"
datadir <- "C:/Users/Rachael/OneDrive - University of Virginia/UVA/Fall_2020/SYS_6021/Project/6021_Project2/"


#sourcedir <- "D:/GoogleDrive/Julie_SYS4021/Projects/TimeSeries2020"
#datadir <- "D:/GoogleDrive/Julie_SYS4021/Data Sets/AirQualityUCI"
opts_knit$set(root.dir = sourcedir)
library(forecast)
library(imputeTS)
library(ggplot2)
library(ggfortify)
library(mtsdi)
library(forecast)
library(lubridate)
library(tidyverse)
library(ggpubr)
library(tseries)
```

# Load data and impute missing values
```{r cars}
setwd(datadir)

airquality = read.csv('AirQualityUCI.csv')

# replace -200 with NA
airquality[airquality == -200] <- NA

# convert integer type to numeric
intcols = c(4,5,7,8,9,10,11,12)
for(i in 1:length(intcols)){
  airquality[,intcols[i]] <- as.numeric(airquality[,intcols[i]])
}

setwd(sourcedir)

# create new data frame with just NO2 and impute missing values
AQdata = airquality["NO2.GT."]
AQdata = na_interpolation(AQdata)

# aggregate to daily maxima for model building
dailyAQ <- aggregate(AQdata, by=list(as.Date(airquality[,1],"%m/%d/%Y")), FUN=max)

# create time series of NO2
NO2.ts <- ts(dailyAQ[,2])

# create time series of NO2 without the last 7 days. (I added this)
dailyAQ.7 <- dailyAQ[1:384,]
NO2.7.ts <- ts(dailyAQ.7[,2])
```

1. Building Univariate Time Series Models: Build a time series model
of daily maximum nitrogen dioxide (NO2) concentrations using all but
the last 7 days of observations. Make sure to address the following items:
(a) How you discovered and modeled any seasonal components, if applicable. (5 points)
(b) How you discovered and modeled any trends, if applicable. (5 points)
(c) How you determined autoregressive and moving average components,
if applicable. (10 points)
(d) How you assessed alternative models (e.g. adjusted R2, AIC, diagnostics, etc.). Assessments should discuss diagnostics and at least one metric. Show and discuss diagnostics of the residuals’ homoscedasticity, Gaussianity, and independence. What problems, if any, remain in the diagnostics of the selected model? (20 points)
(e) Forecast the next 7 days of NO2 concentrations using your selected model. Plot the forecasts vs. true values. What is the MSE of the 7-day forecast? (10 points)

```{r}
# Build a time series model of daily maximum nitrogen dioxide (NO2) concentrations using all but the last 7 days of observations.
# plot the time series
autoplot(NO2.7.ts, ylab="Max NO2 Concentration", xlab="Day")
# seems to be trending upwards and has an increasing variance over time.

oacf <- ggAcf(NO2.7.ts)
opacf <- ggPacf(NO2.7.ts)
oacf
opacf
# acf seems to show a decay over time. Each lag in acf and pacf are significant though and never cut off.
# this tells us that there could be some seasonality maybe?

# Next, model the trend of the concentrations.
time.no2<-c(1:(length(NO2.7.ts)))
time.no27<-c(1:(365))
#Build a new model, spam.trend which predicts spam.ts based on the time variable, time.spam
no27.trend<-lm(NO2.7.ts~time.no2)

##use the summary() command on spam.trend. Is time significant in predicting no2 frequency? yes, it is!
summary(no27.trend)

plot(NO2.7.ts) #linear model
abline(no27.trend,col='red')
# Seems to be an upward trend over time

# Get the periodogram for NO2.7.ts
pg.no27 <- spec.pgram(NO2.7.ts,spans=9,demean=T,log='no')
#spans is ban width, demean is subtract the mean, and log means
# no log transform on the y-axis
#for plotting in ggplot below
spec.no27 <- data.frame(freq=pg.no27$freq, spec=pg.no27$spec)
ggplot(spec.no27) + geom_line(aes(x=freq,y=spec)) + 
  ggtitle("Smooth Periodogram of NO2 Concenctrations")

# There is clearly a peak at the beginning
# Where is the peak?
max.omega.no27<-pg.no27$freq[which(pg.no27$spec==max(pg.no27$spec))]

# Where is the peak?
max.omega.no27
# 0.005208333
# What is the period?
1/max.omega.no27
# period of 192

# Are there any other peaks that we could try and consider in the model?
sorted.spec <- sort(pg.no27$spec, decreasing=T, index.return=T)
names(sorted.spec)

# corresponding periods (omegas = frequences, Ts = periods)
sorted.omegas <- pg.no27$freq[sorted.spec$ix]
sorted.Ts <- 1/pg.no27$freq[sorted.spec$ix]

# look at first 20
sorted.omegas[1:20]
sorted.Ts[1:20]

# Outputs: 192, 128, 384, 96, 76.8, 64, 54.9, 48, 42.7, 38.4
```
```{r}
#*****************************
# 
# Model Trend and Seasonality  
# 
#*****************************
no27.trend<-lm(NO2.7.ts ~ time.no2)

# Use the summary() command on no27.trend,
# Is time significant in predicting the maximum NO2 concentrations
summary(no27.trend)
# concentration coefficeient is not statistically significant, 
#intercept would be the average at time =0
# no, this isn't good enough to see a tredn though
# because we need to account for seasonality

# Plot no27.trend model
ggplot(dailyAQ.7, aes(x=Group.1,y=NO2.GT.)) + geom_line() +
  stat_smooth(method="lm",col="red") + xlab("") + ylab("Richmond Tmin")

# Model diagnostics for no27.trend
autoplot(no27.trend, labels.id = NULL)
    

# Model seasonality, looking at different values from the peak. (# Outputs: 192, 128, 384, 96, 76.8, 64, 54.9, 48, 42.7, 38.4)
no27.trend.seasonal1 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192))
summary(no27.trend.seasonal1)
AIC(no27.trend.seasonal1)
# 1 is significant and positive, AIC of 4000.19

# Try including another term, 128
no27.trend.seasonal2 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128))
summary(no27.trend.seasonal2)
AIC(no27.trend.seasonal2)
# All trig terms are significant, AIC of 3984.142. Therefore, this one is better than the prior.

# Try including a third, 384
no27.trend.seasonal3 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/384) + cos(2*pi*time.no2/384))
summary(no27.trend.seasonal3)
AIC(no27.trend.seasonal3)
# All trig terms are significant again, AIC of 3912.906. Therefore the 2nd is still the best.

# Try a different third term, 96.
no27.trend.seasonal4 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/96) + cos(2*pi*time.no2/96))
summary(no27.trend.seasonal4)
AIC(no27.trend.seasonal4)
# The 96 trig terms are not significant, and AIC of 3985.905, so it's closer than the prior. 

# Try another from the list as a 3rd, 76.8
no27.trend.seasonal5 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8))
summary(no27.trend.seasonal5)
AIC(no27.trend.seasonal5)
# Each one has 1 significant trig term, AIC of 3961.554. This is now our best.

# Try adding a fourth term, 64
no27.trend.seasonal6 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64))
summary(no27.trend.seasonal6)
AIC(no27.trend.seasonal6)
# Each has at least 1 significant trig term, AIC of 3925.592, our new best.

# Try adding a fifth term, 54.9
no27.trend.seasonal7 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64)+sin(2*pi*time.no2/54.9) + cos(2*pi*time.no2/54.9))
summary(no27.trend.seasonal7)
AIC(no27.trend.seasonal7)
# Not every term has at least 1 significant trig.

# Try a different fifth term, 48
no27.trend.seasonal8 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64)+sin(2*pi*time.no2/48) + cos(2*pi*time.no2/48))
summary(no27.trend.seasonal8)
AIC(no27.trend.seasonal8)
# Not every term has at least 1 significant trig.

# Try a different fifth term, 42.7
no27.trend.seasonal9 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64)+sin(2*pi*time.no2/42.7) + cos(2*pi*time.no2/42.7))
summary(no27.trend.seasonal9)
AIC(no27.trend.seasonal9)
# Not every term has at least 1 significant trig.

# Try a different fifth term, 38.4
no27.trend.seasonal10 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64)+sin(2*pi*time.no2/38.4) + cos(2*pi*time.no2/38.4))
summary(no27.trend.seasonal10)
AIC(no27.trend.seasonal10)
# Not every term has at least 1 significant trig. Therefore, the 6th is still the best. 

# We'll try another variation of the 6th with possible patterns we can think of, how about 30, 7, and 365 (months, weeks, and years)
no27.trend.seasonal6.2 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64)+ sin(2*pi*time.no2/7) + cos(2*pi*time.no2/7))
summary(no27.trend.seasonal6.2)
AIC(no27.trend.seasonal6.2)
# With 7, all have at least 1 significant trig term, and the lowest AIC pf 3895.892

# Try adding 30 to the above.
no27.trend.seasonal6.2.2 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64)+ sin(2*pi*time.no2/7) + cos(2*pi*time.no2/7)+sin(2*pi*time.no2/30) + cos(2*pi*time.no2/30))
summary(no27.trend.seasonal6.2.2)
AIC(no27.trend.seasonal6.2.2)
# With 30, all have at least 1 significant trig term, with the best AIC of 3889.858.

# Try adding 365 to the above.
no27.trend.seasonal6.2.2.2 <- lm(NO2.7.ts ~ time.no2 + sin(2*pi*time.no2/192) + cos(2*pi*time.no2/192) + sin(2*pi*time.no2/128) + cos(2*pi*time.no2/128)+ sin(2*pi*time.no2/76.8) + cos(2*pi*time.no2/76.8)+ sin(2*pi*time.no2/64) + cos(2*pi*time.no2/64)+ sin(2*pi*time.no2/7) + cos(2*pi*time.no2/7)+sin(2*pi*time.no2/30) + cos(2*pi*time.no2/30)+sin(2*pi*time.no2/365) + cos(2*pi*time.no2/365))
summary(no27.trend.seasonal6.2.2.2)
AIC(no27.trend.seasonal6.2.2.2)
# With 30, all have at least 1 significant trig term, with the best AIC of 3804.863
# Let's use model 6.2.2.2, we will rename it to 
no27.trend.seasonal <-
  no27.trend.seasonal6.2.2.2
# now we know concentration is significant and it is positive. Note the cos term is not significant though.
# Plot no27.trend.seasonal model
ggplot(dailyAQ.7, aes(x=Group.1,y=NO2.GT.)) + geom_line() + 
  geom_line(aes(x=Group.1,y=no27.trend.seasonal$fitted.values),color="red") 

# Get the residuals from the no27.trend.seasonal model above and store in e.ts.no27:
e.ts.no27 <- ts(no27.trend.seasonal$residuals)
    
# Plot the residuals for the no27.trend model
autoplot(e.ts.no27)

# we want it to have an average mean of constant 0
# it is close, but there still is some non-independence
# that we need to capture as well. There seems to be some inconsistent extremes

# Plot the autocorrelation (ACF) of the residuals of no27.trend.seasonal
no27.acf <- ggAcf(e.ts.no27)
no27.acf
# several lags are significant, there is some sinosoidal
# behavior, look at PACF next

# Plot the partial autocorrelation (PACF) of the residuals no27.trend.seasonal
no27.pacf <- ggPacf(e.ts.no27)
no27.pacf     
# also has sinusoidal behavior and doesn't cut off
# which suggests we may have autoregressive and moving average terms

# We could have some autoregressive terms to be examined.
# Choose p and q terms for e.ts.no27 based on the acf and pacf. Doesn't seem to have linear decay so no need to check differencing terms but will check a first order differencing to be sure.
no27.diff.acf <- ggAcf(diff(e.ts.no27))
no27.diff.pacf <- ggPacf(diff(e.ts.no27))
ggarrange(no27.diff.acf,no27.diff.pacf,nrow=2,ncol=1)
# These appear to be worse, so we were right.

# therefore, the fact that ARMA has best aic seems to show that
# our hypothesis that it has autoregressive and moving average terms
#use automated selection method to see what p and q values are the best

# using arima() function to find best p, q, and d values
no27.auto <- auto.arima(e.ts.no27, approximation= FALSE)
summary(no27.auto)
#(1, 0, 0) (p, d, q), 
# aic = 3707.15
```
```{r}
# Examining based off of TS3 lecture:

# ar(2) p=2 (try 2 because auto has 1)
no27.ar1 <- arima(e.ts.no27, order=c(2,0,0), include.mean=FALSE)
summary(no27.ar1)     #3709.15

# ma(2) p=0, q=2
no27.ma2 <- arima(e.ts.no27, order=c(0,0,2), include.mean=FALSE)
summary(no27.ma2)     #3712.34

# arma(1,2) p=1, q=2
no27.arma12 <- arima(e.ts.no27, order=c(1,0,2), include.mean=FALSE)
summary(no27.arma12)     #3710.73

# auto from above
summary(no27.auto)     #3707.15
# Based on AIC, no27.auto appears to be the best

# BIC (penalizes additional terms more)
BIC(no27.ar1)     #3721.001
BIC(no27.ma2)     #3724.187
BIC(no27.arma12)     #3726.529
BIC(no27.auto)     #3715.053
# Based on BIC. no27.auto is still the best 

```
```{r}
# assess residuals vs. fitted
model1 = ggplot() + geom_point(aes(x=fitted(no27.ar1), y=no27.ar1$residuals)) + ggtitle("AR1")
model2 = ggplot() + geom_point(aes(x=fitted(no27.ma2), y=no27.ma2$residuals)) + ggtitle("MA2")
model3 = ggplot() + geom_point(aes(x=fitted(no27.arma12), y=no27.arma12$residuals)) + ggtitle("ARMA12")
model4 = ggplot() + geom_point(aes(x=fitted(no27.auto), y=no27.auto$residuals)) + ggtitle("Auto")

ggarrange(model1, model2, model3, model4, ncol=2, nrow=2)
# all seem to have constant mean of zero, no significant differences


# assess normality of residuals
model1 = qplot(sample=no27.ar1$residuals) + stat_qq_line(color="red") + ggtitle("AR1")
model2 = qplot(sample=no27.ma2$residuals) + stat_qq_line(color="red") + ggtitle("MA2")
model3 = qplot(sample=no27.arma12$residuals) + stat_qq_line(color="red") + ggtitle("ARMA12")
model4 = qplot(sample=no27.auto$residuals) + stat_qq_line(color="red") + ggtitle("Auto")

ggarrange(model1, model2, model3, model4, ncol=2, nrow=2)
# all look roughly the same and pretty good, looks pretty gaussian
# so we're now mostly interested in seeing if our residuals are independent

# Plot diagnostics for independence of residuals using tsdiag()
ggtsdiag(no27.ar1,gof.lag=20)
# good for up to 5 lags
ggtsdiag(no27.ma2,gof.lag=20)
# up to 2 lags
ggtsdiag(no27.arma12,gof.lag=20)
# about 6 lags
ggtsdiag(no27.auto,gof.lag=20)
# good for 4 lags

# no27arma12, no27.auto, and no27.ar1 are the most adequate, while
# no27.ma2 isn't great

# Plot the autocorrelation (ACF) and partial autocorrelation (PACF) of the residuals of no27.auto
no27.auto.resid.acf <- ggAcf(no27.auto$residuals)
no27.auto.resid.pacf <- ggPacf(no27.auto$residuals)
ggarrange(no27.auto.resid.acf, no27.auto.resid.pacf,nrow=2,ncol=1)

# The no27.auto model accounts for correlation in the residuals, 5

# Plot the autocorrelation (ACF) and partial autocorrelation (PACF) of the residuals of no27.arma12
no27.arma12.resid.acf <- ggAcf(no27.arma12$residuals)
no27.arma12.resid.pacf <- ggPacf(no27.arma12$residuals)
ggarrange(no27.arma12.resid.acf, no27.arma12.resid.pacf,nrow=2,ncol=1)
# Does account for the correlation in the residuals, about 5

# Plot the autocorrelation (ACF) and partial autocorrelation (PACF) of the residuals of no27.arma12
no27.ar1.resid.acf <- ggAcf(no27.ar1$residuals)
no27.ar1.resid.pacf <- ggPacf(no27.ar1$residuals)
ggarrange(no27.ar1.resid.acf, no27.ar1.resid.pacf,nrow=2,ncol=1)
# Does account for the correlation in the residuals, about 5

#3 Double check, how to tell if accounts for residuals 

# Plot the best model's fitted values vs. true values
Tmin.fitted <- no27.trend.seasonal$fitted.values + fitted(no27.auto)

ggplot() + geom_line(aes(x=time.no2,y=NO2.7.ts[1:length(time.no2)],color="True")) +
  geom_line(aes(x=time.no2,y=Tmin.fitted,color="Fitted")) + xlab("Time") + 
  ylab("Max NO2 Concentrations")
# Looks pretty good, some of the peaks are off, but the patterns are similar.
```

```{r}
# Forecasting the next 7 days with the best model: (using auto for now)

no27.auto.forecast <- forecast(no27.auto, h=7)
autoplot(no27.auto.forecast,main="Forecasts from AUTO(1,0,0) with zero mean")

# The actual time series for the period including the last week is NO2.ts
# Prediction for the next week by no27.auto:
next.wk.time <- c((length(NO2.ts)-6):(length(NO2.ts)))

next.wk <- data.frame(time.no2 = next.wk.time, no27 = NO2.ts[next.wk.time])

# The actual time series for the test period
next.wk.ts <- NO2.ts[next.wk.time]
next.wk.ts <- ts(next.wk$no27)

E_Y.pred <- predict(no27.trend.seasonal, newdata=next.wk)
e_t.pred <- forecast(no27.auto, h=7)
next.wk.prediction <- E_Y.pred + e_t.pred$mean

# MSE:

mean((next.wk.prediction - next.wk$no27)^2)     #345.7277

# Plot actual values and predicted values
plot(ts(next.wk$no27),type='o',ylim=c(110,300))
lines(ts(next.wk.prediction),col='red',type='o')
lines(1:7, E_Y.pred + e_t.pred$lower[,2], col = "red", lty = "dashed")
lines(1:7, E_Y.pred + e_t.pred$upper[,2], col = "red", lty = "dashed")
#legend(1,200, legend = c("Actual", "Predicted"), lwd = 2, col = c("black", "red")) 

# Note: Can repeat this process with other models as needed.
```

# End of part 1!

2. Simulating Univariate Time Series Models: Simulate a year of synthetic observations of daily maximum nitrogen dioxide (NO2) concentrations from your selected model. Set the seed so you will get the same
results each time. You will need to consider the sum of the linear models
of the trend + seasonality, and the residual models. Assess and compare
the model’s performance with respect to: (50 points)

(a) Ability to reproduce appearance of time series. Plot observations and
simulations and visually compare their characteristics. (10 points)

(b) Ability to reproduce observed trends. You can assess this by building
a linear model of the trend + seasonality of the simulations and comparing the coefficient estimates with the linear model of the trend +
seasonality of the observations. What is the percent difference in the
coefficient on time? (10 points)

(c) Ability to reproduce seasonality of the time series. Analysis can be
visual, simply comparing the periodogram of the observations and
simulations. (10 points)

(d) Ability to reproduce observed mean and variance of the time series
(Hint: Use the functions ‘mean(ts)’ and ‘var(ts)’ where ts is a time
series, and find the percent difference between observations and simulations) (10 points)

(e) Ability to reproduce the autocorrelation of the time series. Analysis
can be visual, simply comparing the ACF and PACF of the observations and simulations. (10 points)

```{r}
# Simulate a year of daily maximum NO2 concentrations with the best model
set.seed(1)
auto.sim <- arima.sim(n=365, list(ar=c(no27.auto$coef[1])), sd=sqrt(no27.auto$sigma2))
# Add mean predictions and plot simulation of Tmax concentrations
next.yr.time <- c(1:(365))
next.yr <- data.frame(time.no2 = next.yr.time)

next.yr.predictions <- predict(no27.trend.seasonal, newdata = next.yr)

# plot simulated concentrations
autoplot(ts(next.yr.predictions + auto.sim),xlab="Time",ylab="Simulated Daily Maximum Concentrations")

# Now plot simulated vs actual
dailyAQ.365 <- dailyAQ[1:365,]
NO2.365.ts <- ts(dailyAQ.365[,2])
sim =autoplot(ts(next.yr.predictions + auto.sim),xlab="Time",ylab="Simulated Daily Maximum Concentrations")
act =autoplot(NO2.365.ts, xlab="Time", ylab="Actual NO2 level")
act
sim
ggarrange(sim, act ,nrow=2,ncol=1)

# Next, look at periodogram
sim.ts <- next.yr.predictions + auto.sim
as <- pg.no27 #actual periodogram formed earlier
ps <- spec.pgram(sim.ts,spans=9,demean=T,log='no')

spec.no27 <- data.frame(freq=pg.no27$freq, spec=pg.no27$spec)
ggplot(spec.no27) + geom_line(aes(x=freq,y=spec)) + 
  ggtitle("Smooth Periodogram of NO2 Concenctrations")
pspec.no27 <- data.frame(freq=ps$freq, spec=ps$spec)
ggplot(pspec.no27) + geom_line(aes(x=freq,y=spec)) + 
  ggtitle("Smooth Periodogram of Predicted NO2 Concenctrations")
```

```{r}
# Model the simulated data and see how the trend, seasonality, and residuals look.
sim <- as.numeric(sim.ts)

simno27.ts <- lm(sim~next.yr.time+sin(2*pi*next.yr.time/192)
                    +cos(2*pi*next.yr.time/192)+sin(2*pi*next.yr.time/128)
                    +cos(2*pi*next.yr.time/128)+sin(2*pi*next.yr.time/76.8)
                    +cos(2*pi*next.yr.time/76.8)+sin(2*pi*next.yr.time/64)
                    +cos(2*pi*next.yr.time/64)+sin(2*pi*next.yr.time/64)
                    +cos(2*pi*next.yr.time/64)+sin(2*pi*next.yr.time/7)
                    +cos(2*pi*next.yr.time/7)+sin(2*pi*next.yr.time/30)
                    +cos(2*pi*next.yr.time/30)+sin(2*pi*next.yr.time/365)
                    +cos(2*pi*next.yr.time/365))

summary(simno27.ts)     # Each has at least one significant trig term again.
AIC(simno27.ts)     #3556.824

no27.trend.seasonal$coefficients[2]
simno27.ts$coefficients[2]

# Look at the differences in the first coefficient (ignore intercept here)
abs(no27.trend.seasonal$coefficients[2]-simno27.ts$coefficients[2])/no27.trend.seasonal$coefficients[2]     #1.018155 

ggplot(sim.ts, aes(x=next.yr.time,y=sim)) + geom_line() +   geom_line(aes(x=next.yr.time,y=simno27.ts$fitted.values),color="red") + xlab("") + ylab("NO2")

```

```{r}
# Use the functions ‘mean(ts)’ and ‘var(ts)’ where ts is a time series, and find the percent difference between observations and simulations.

# Means
a <- mean(sim.ts)
b <- mean(NO2.7.ts)
abs(b-a)/b     #0.00423076

# Variances
c <- var(sim.ts)
d <- var(NO2.7.ts)
abs(d-c)/d     #0.1139395

```
```{r}
# Compare the ACF and PACF of the observations and simulations
sim.res <- simno27.ts$residuals
sim.res.ts <- ts(sim.res)

sim.res.acf <- ggAcf(sim.res.ts)
sim.res.pacf <- ggPacf(sim.res.ts)

ggarrange(oacf, sim.res.acf, opacf, sim.res.pacf, nrow=2, ncol=2)

# ACF looks very different, pacf looks more similar
```