.Rhistory

y = "Residuals")
poly_season_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_season)) +
geom_point() +
labs(title = "Season Dummy Variable Third Order Polynomial Residual",
# subtitle = "Model includes season as a categorical variable (Winter, Spring, Summer, Autumn)",
x = "Year",
y = "Residuals")
# Combine the three plots side by side
residual_combined_plot <- poly_month_plot + poly_month_year_plot + poly_season_plot +
plot_layout(nrow = 3) # Set layout with 3 columns
# Display the combined plot
residual_combined_plot
rm(list = ls())
library(tidyverse)
library(tsibble)
library(latex2exp)
library(feasts)
library(patchwork)
library(forecast)
library(tseries)
library(knitr)
theme_set(theme_bw())
knitr::opts_chunk$set(dpi=1000)
co2_tsib<-as_tsibble(co2)
co2_tsib %>%
ggplot() +
aes(x=index, y=value) +
geom_line(color = 'steelblue') +
labs(
title = TeX(r'(Monthly Mean $CO_2$)'),
subtitle = 'The "Keeling Curve"',
x = 'Month and Year',
y = TeX(r'($CO_2$ parts per million)')
)
co2_tsib
# monthly time series with the line of best fit
co2_trend_plot <- co2_tsib %>%
ggplot(aes(x = index, y = value)) +
geom_line(color = 'black', size = .5) +
geom_smooth(method = "lm", formula = "y ~ x",se = F, color = 'blue', size = .8) +
labs(title = 'Average Trend Plotted on Monthly Average Co2 ppmv') +
xlab('Year') +
ylab('Co2 concentrations (ppmv)')
# average_yearly_increase
co2_tsib_yearly_change <- co2_tsib %>% as_tibble() %>%
mutate(year = year(index)) %>%
group_by(year) %>%
summarise(`yearly_co2` = mean(value)) %>%
ungroup() %>%
mutate(lag_co2 = lag(yearly_co2),
change = yearly_co2 - lag_co2,
percent_change = ((yearly_co2 - lag_co2)/yearly_co2)*100)
# getting average increase (i.e. size of the trend)
yearly_mean <- mean(co2_tsib_yearly_change$change, na.rm = T) # average 1.26 units of co2 change each year
yearly_sd <- sd(co2_tsib_yearly_change$change, na.rm = T) # with a sd of .51 units of co2
change_hist <- ggplot(co2_tsib_yearly_change, aes(x = change)) +
geom_histogram(color = 'gray20', fill = 'gray', binwidth = .25) +
scale_x_continuous(breaks = seq(floor(min(co2_tsib_yearly_change$change, na.rm =T)),
ceiling(max(co2_tsib_yearly_change$change, na.rm =T)), by = 0.5)) +
ggtitle('Histogram of Yearly changes in Co2 ppmv')
# bulk of increases yearly seem to be between .5 and 2 units on the
co2_trend_plot / change_hist
# inspecting acf and graph of co2 concentrations over time
co2_acf <- acf(co2_tsib$value, plot = F)
co2_acf_plot <-  autoplot(co2_acf) +
labs(title = "ACF plot of monthly co2 Concentrations", x = 'lag', y = 'Autocorrelation')
monthly_co2_ave_plot <- co2_tsib %>% as_tibble() %>%
mutate(month = month(index)) %>%
group_by(month) %>%
summarise(co2_monthly_ave = mean(value, na.rm = T)) %>%
mutate(month_str = factor(month.abb[month], levels = month.abb)) %>%
ungroup() %>%
ggplot(aes(x = month_str, y = co2_monthly_ave, group = 1)) +
geom_line(size = .8, color = 'purple4') +
geom_point(size = 1.5) +
ggtitle("Average Co2 Concentrations Across\nEach Month") +
xlab('Month') +
ylab('Co2 Concentrations (ppmv)') +
theme(axis.text.x = element_text(angle = 45))
co2_acf_plot | monthly_co2_ave_plot
# making a plot to show how the relationship looks like with yearly averages over the seasons
yearly_ave_w_residuals <- co2_tsib %>% as_tibble() %>%
mutate(year = year(index)) %>%
group_by(year) %>%
mutate(`Yearly Co2` = mean(value)) %>%
mutate(residual = value - `Yearly Co2`) %>%
pivot_longer(cols = c(value, `Yearly Co2`, residual), names_to = "type", values_to = "Monthly Co2") %>%
mutate(residual_bool = if_else(type == "residual", "Residuals", "Monthly Time Series Plotted on Yearly Average Co2"))
yearly_ave_w_residuals_plot <- yearly_ave_w_residuals %>%
ggplot(aes(x = index, y = `Monthly Co2`, color = type)) +
geom_line() +
facet_wrap(~residual_bool, scales = "free_y", ncol = 1) +
xlab('Date') +
theme(legend.position = "none")
# residuals of a simple yearly average look fairly stationary
# with some years having larger and smaller co2 variances
adf_result <- yearly_ave_w_residuals %>%
filter(residual_bool == 'Residuals') %>%
ungroup() %>%
pull(`Monthly Co2`) %>%
adf.test()
# Extract the relevant results into a data frame
adf_summary <- data.frame(
Statistic = adf_result$statistic,
P_Value = adf_result$p.value,
Method = adf_result$method,
Alternative = adf_result$alternative
)
yearly_ave_w_residuals_plot
kable(adf_summary, caption = "ADF Test Results")
# monthly time series with the curvilinear line of best fit
co2_curv_trend_plot <-
co2_tsib %>%
ggplot(aes(x = index, y = value)) +
geom_line(color = 'black', size = .5) +
geom_smooth(method = "lm", formula = "y ~ poly(x,3)",se = F, color = 'blue', size = .8) +
labs(title = 'Curvilinear Trend Plotted on Monthly Average Co2 ppmv') +
xlab('Year') +
ylab('Co2 concentrations (ppmv)')
co2_trend_plot / co2_curv_trend_plot
co2_linear_model <- lm(value ~ index, data = co2_tsib)
# Plot residuals
co2_tsib$residuals_linear <- residuals(co2_linear_model)
lm_residuals_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_linear)) +
geom_point() +
labs(title = "Residuals of Linear Time Trend Model", x = "Year", y = "Residuals")
lm_residuals_plot
co2_quad_model <- lm(value ~ poly(index, 2), data = co2_tsib)
# Plot residuals
co2_tsib$residuals_quad <- residuals(co2_quad_model)
quad_residuals_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_quad)) +
geom_point() +
labs(title = "Residuals of Quadratic Time Trend Model", x = "Year", y = "Residuals")
quad_residuals_plot
co2_poly_model <- lm(value ~ poly(index,3), data = co2_tsib)
# Plot residuals
co2_tsib$residuals_poly <- residuals(co2_poly_model)
poly_residuals_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly)) +
geom_point() +
labs(title = "Residuals of Polynomial Time Trend Model", x = "Year", y = "Residuals")
poly_residuals_plot
# Log transformation
co2_tsib$log_value <- log(co2_tsib$value)
# Plot log-transformed data
log_data_plot <- ggplot(co2_tsib, aes(x = index, y = log_value)) +
geom_line() +
labs(title = "Log-transformed Co2 Concentrations", x = "Year", y = "Log of Co2 concentrations")
log_data_plot
# Create seasonal dummy variables
co2_tsib$month <- factor(month(co2_tsib$index))
co2_tsib$year <- factor(year(co2_tsib$index))
# Define a function to convert months into seasons
co2_tsib$season <- case_when(
month(co2_tsib$index) %in% c(12, 1, 2) ~ "Winter",
month(co2_tsib$index) %in% c(3, 4, 5) ~ "Spring",
month(co2_tsib$index) %in% c(6, 7, 8) ~ "Summer",
month(co2_tsib$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Convert season into a factor
co2_tsib$season <- factor(co2_tsib$season, levels = c("Winter", "Spring", "Summer", "Autumn"))
# Fit polynomial model with seasonal dummies
poly_month <- lm(value ~ poly(index, 3) + month, data = co2_tsib)
poly_month_year <- lm(value ~ poly(index, 3) + month + year, data = co2_tsib)
poly_season <- lm(value ~ poly(index, 3) + season, data = co2_tsib)
# Summary of the model
summary(poly_seasonal_model)
# Create seasonal dummy variables
co2_tsib$month <- factor(month(co2_tsib$index))
co2_tsib$year <- factor(year(co2_tsib$index))
# Define a function to convert months into seasons
co2_tsib$season <- case_when(
month(co2_tsib$index) %in% c(12, 1, 2) ~ "Winter",
month(co2_tsib$index) %in% c(3, 4, 5) ~ "Spring",
month(co2_tsib$index) %in% c(6, 7, 8) ~ "Summer",
month(co2_tsib$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Convert season into a factor
co2_tsib$season <- factor(co2_tsib$season, levels = c("Winter", "Spring", "Summer", "Autumn"))
# Fit polynomial model with seasonal dummies
poly_month <- lm(value ~ poly(index, 3) + month, data = co2_tsib)
poly_month_year <- lm(value ~ poly(index, 3) + month + year, data = co2_tsib)
poly_season <- lm(value ~ poly(index, 3) + season, data = co2_tsib)
# Get residuals
co2_tsib$residuals_poly_month <- residuals(poly_month)
co2_tsib$residuals_poly_month_year <- residuals(poly_month_year)
co2_tsib$residuals_poly_season <- residuals(poly_season)
# residual plots
poly_month_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month)) +
geom_point() +
labs(title = "Month Dummy Variable Third Order Polynomial Residuals",
# subtitle = "Model includes month as a categorical variable",
x = "Year",
y = "Residuals")
poly_month_year_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month_year)) +
geom_point() +
labs(title = "Month and Year Dummy Variables Third Order Polynomial Residual",
# subtitle = "Model includes both month and year as categorical variables",
x = "Year",
y = "Residuals")
poly_season_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_season)) +
geom_point() +
labs(title = "Season Dummy Variable Third Order Polynomial Residual",
# subtitle = "Model includes season as a categorical variable (Winter, Spring, Summer, Autumn)",
x = "Year",
y = "Residuals")
# Combine the three plots side by side
residual_combined_plot <- poly_month_plot + poly_month_year_plot + poly_season_plot +
plot_layout(nrow = 3) # Set layout with 3 columns
# Display the combined plot
print(residual_combined_plot)
summary(poly_month)
summary(poly_month_year)
summary(poly_season)
# Create seasonal dummy variables
co2_tsib$month <- factor(month(co2_tsib$index))
co2_tsib$year <- factor(year(co2_tsib$index))
# Define a function to convert months into seasons
co2_tsib$season <- case_when(
month(co2_tsib$index) %in% c(12, 1, 2) ~ "Winter",
month(co2_tsib$index) %in% c(3, 4, 5) ~ "Spring",
month(co2_tsib$index) %in% c(6, 7, 8) ~ "Summer",
month(co2_tsib$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Convert season into a factor
co2_tsib$season <- factor(co2_tsib$season, levels = c("Winter", "Spring", "Summer", "Autumn"))
# Fit polynomial model with seasonal dummies
poly_month <- lm(value ~ poly(index, 3) + month, data = co2_tsib)
poly_month_year <- lm(value ~ poly(index, 3) + month + year, data = co2_tsib)
poly_season <- lm(value ~ poly(index, 3) + season, data = co2_tsib)
# Get residuals
co2_tsib$residuals_poly_month <- residuals(poly_month)
co2_tsib$residuals_poly_month_year <- residuals(poly_month_year)
co2_tsib$residuals_poly_season <- residuals(poly_season)
# residual plots
poly_month_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month)) +
geom_point() +
labs(title = "Month Dummy Variable Third Order Polynomial Residuals",
# subtitle = "Model includes month as a categorical variable",
x = "Year",
y = "Residuals")
poly_month_year_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month_year)) +
geom_point() +
labs(title = "Month and Year Dummy Variables Third Order Polynomial Residual",
# subtitle = "Model includes both month and year as categorical variables",
x = "Year",
y = "Residuals")
poly_season_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_season)) +
geom_point() +
labs(title = "Season Dummy Variable Third Order Polynomial Residual",
# subtitle = "Model includes season as a categorical variable (Winter, Spring, Summer, Autumn)",
x = "Year",
y = "Residuals")
# Combine the three plots side by side
residual_combined_plot <- poly_month_plot + poly_month_year_plot + poly_season_plot +
plot_layout(nrow = 3) # Set layout with 3 columns
# Display the combined plot
print(residual_combined_plot)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index) + 1/12, to = 2020, by = 1/12)
# Ensure index is in Date format
co2_tsib$index_date <- as.Date(co2_tsib$index)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + 1/12, to = 2020, by = 1/12)
# Ensure index is in Date format
co2_tsib$index_date <- as.Date(co2_tsib$index)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + months(1), to = as.Date("2020-12-01"), by = "month")
future_data <- data.frame(index = future_years, month = factor(month(future_years)))
# Predict using the model
future_data$forecast <- predict(poly_season, newdata = future_data)
# Ensure index is in Date format
co2_tsib$index_date <- as.Date(co2_tsib$index)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + months(1), to = as.Date("2020-12-01"), by = "month")
future_data <- data.frame(index = future_years, month = factor(month(future_years)))
# Get the season
future_data$season <- case_when(
month(future_data$index) %in% c(12, 1, 2) ~ "Winter",
month(future_data$index) %in% c(3, 4, 5) ~ "Spring",
month(future_data$index) %in% c(6, 7, 8) ~ "Summer",
month(future_data$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Predict using the model
future_data$forecast <- predict(poly_season, newdata = future_data)
# Plot the forecasts
forecast_plot <- ggplot() +
geom_line(data = co2_tsib, aes(x = index, y = value), color = "black") +
geom_line(data = future_data, aes(x = index, y = forecast), color = "red") +
labs(title = "Forecast of Co2 Concentrations to 2020", x = "Year", y = "Co2 concentrations (ppmv)")
# Ensure index is in Date format
co2_tsib$index_date <- as.Date(co2_tsib$index)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + months(1), to = as.Date("2020-12-01"), by = "month")
future_data <- data.frame(index = future_years, month = factor(month(future_years)))
# Get the season
future_data$season <- case_when(
month(future_data$index) %in% c(12, 1, 2) ~ "Winter",
month(future_data$index) %in% c(3, 4, 5) ~ "Spring",
month(future_data$index) %in% c(6, 7, 8) ~ "Summer",
month(future_data$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Predict using the model
future_data$forecast <- predict(poly_season, newdata = future_data)
# Plot the forecasts
forecast_plot <- ggplot() +
geom_line(data = co2_tsib, aes(x = index, y = value), color = "black") +
geom_line(data = future_data, aes(x = index, y = forecast), color = "red") +
labs(title = "Forecast of Co2 Concentrations to 2020", x = "Year", y = "Co2 concentrations (ppmv)")
forecast_plot
View(future_data)
# Ensure index is in Date format
co2_tsib$index_date <- as.Date(co2_tsib$index)
# Create seasonal dummy variables
co2_tsib$month <- factor(month(co2_tsib$index))
co2_tsib$year <- factor(year(co2_tsib$index))
# Define a function to convert months into seasons
co2_tsib$season <- case_when(
month(co2_tsib$index) %in% c(12, 1, 2) ~ "Winter",
month(co2_tsib$index) %in% c(3, 4, 5) ~ "Spring",
month(co2_tsib$index) %in% c(6, 7, 8) ~ "Summer",
month(co2_tsib$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Convert season into a factor
co2_tsib$season <- factor(co2_tsib$season, levels = c("Winter", "Spring", "Summer", "Autumn"))
# Fit polynomial model with seasonal dummies
poly_month <- lm(value ~ poly(index_date, 3) + month, data = co2_tsib)
poly_month_year <- lm(value ~ poly(index_date, 3) + month + year, data = co2_tsib)
poly_season <- lm(value ~ poly(index_date, 3) + season, data = co2_tsib)
# Get residuals
co2_tsib$residuals_poly_month <- residuals(poly_month)
co2_tsib$residuals_poly_month_year <- residuals(poly_month_year)
co2_tsib$residuals_poly_season <- residuals(poly_season)
# residual plots
poly_month_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month)) +
geom_point() +
labs(title = "Month Dummy Variable Third Order Polynomial Residuals",
# subtitle = "Model includes month as a categorical variable",
x = "Year",
y = "Residuals")
poly_month_year_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month_year)) +
geom_point() +
labs(title = "Month and Year Dummy Variables Third Order Polynomial Residual",
# subtitle = "Model includes both month and year as categorical variables",
x = "Year",
y = "Residuals")
poly_season_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_season)) +
geom_point() +
labs(title = "Season Dummy Variable Third Order Polynomial Residual",
# subtitle = "Model includes season as a categorical variable (Winter, Spring, Summer, Autumn)",
x = "Year",
y = "Residuals")
# Combine the three plots side by side
residual_combined_plot <- poly_month_plot + poly_month_year_plot + poly_season_plot +
plot_layout(nrow = 3) # Set layout with 3 columns
# Display the combined plot
print(residual_combined_plot)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + months(1), to = as.Date("2020-12-01"), by = "month")
future_data <- data.frame(index = future_years, month = factor(month(future_years)))
# Get the season
future_data$season <- case_when(
month(future_data$index) %in% c(12, 1, 2) ~ "Winter",
month(future_data$index) %in% c(3, 4, 5) ~ "Spring",
month(future_data$index) %in% c(6, 7, 8) ~ "Summer",
month(future_data$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Predict using the model
future_data$forecast <- predict(poly_month, newdata = future_data)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + months(1), to = as.Date("2020-12-01"), by = "month")
future_data <- data.frame(index = future_years, month = factor(month(future_years)))
# Create seasonal dummy variables
future_data$month <- factor(month(future_data$index))
# Get the season
future_data$season <- case_when(
month(future_data$index) %in% c(12, 1, 2) ~ "Winter",
month(future_data$index) %in% c(3, 4, 5) ~ "Spring",
month(future_data$index) %in% c(6, 7, 8) ~ "Summer",
month(future_data$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Predict using the model
future_data$forecast <- predict(poly_month, newdata = future_data)
# Ensure index is in Date format
co2_tsib$index_date <- as.Date(co2_tsib$index)
# Create seasonal dummy variables
co2_tsib$month <- factor(month(co2_tsib$index))
co2_tsib$year <- factor(year(co2_tsib$index))
# Define a function to convert months into seasons
co2_tsib$season <- case_when(
month(co2_tsib$index) %in% c(12, 1, 2) ~ "Winter",
month(co2_tsib$index) %in% c(3, 4, 5) ~ "Spring",
month(co2_tsib$index) %in% c(6, 7, 8) ~ "Summer",
month(co2_tsib$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Convert season into a factor
co2_tsib$season <- factor(co2_tsib$season, levels = c("Winter", "Spring", "Summer", "Autumn"))
# Fit polynomial model with seasonal dummies
poly_month <- lm(value ~ poly(index 3) + month, data = co2_tsib)
# Ensure index is in Date format
co2_tsib$index_date <- as.Date(co2_tsib$index)
# Create seasonal dummy variables
co2_tsib$month <- factor(month(co2_tsib$index))
co2_tsib$year <- factor(year(co2_tsib$index))
# Define a function to convert months into seasons
co2_tsib$season <- case_when(
month(co2_tsib$index) %in% c(12, 1, 2) ~ "Winter",
month(co2_tsib$index) %in% c(3, 4, 5) ~ "Spring",
month(co2_tsib$index) %in% c(6, 7, 8) ~ "Summer",
month(co2_tsib$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Convert season into a factor
co2_tsib$season <- factor(co2_tsib$season, levels = c("Winter", "Spring", "Summer", "Autumn"))
# Fit polynomial model with seasonal dummies
poly_month <- lm(value ~ poly(index, 3) + month, data = co2_tsib)
poly_month_year <- lm(value ~ poly(index, 3) + month + year, data = co2_tsib)
poly_season <- lm(value ~ poly(index, 3) + season, data = co2_tsib)
# Get residuals
co2_tsib$residuals_poly_month <- residuals(poly_month)
co2_tsib$residuals_poly_month_year <- residuals(poly_month_year)
co2_tsib$residuals_poly_season <- residuals(poly_season)
# residual plots
poly_month_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month)) +
geom_point() +
labs(title = "Month Dummy Variable Third Order Polynomial Residuals",
# subtitle = "Model includes month as a categorical variable",
x = "Year",
y = "Residuals")
poly_month_year_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_month_year)) +
geom_point() +
labs(title = "Month and Year Dummy Variables Third Order Polynomial Residual",
# subtitle = "Model includes both month and year as categorical variables",
x = "Year",
y = "Residuals")
poly_season_plot <- ggplot(co2_tsib, aes(x = index, y = residuals_poly_season)) +
geom_point() +
labs(title = "Season Dummy Variable Third Order Polynomial Residual",
# subtitle = "Model includes season as a categorical variable (Winter, Spring, Summer, Autumn)",
x = "Year",
y = "Residuals")
# Combine the three plots side by side
residual_combined_plot <- poly_month_plot + poly_month_year_plot + poly_season_plot +
plot_layout(nrow = 3) # Set layout with 3 columns
# Display the combined plot
print(residual_combined_plot)
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + months(1), to = as.Date("2020-12-01"), by = "month")
future_data <- data.frame(index = future_years, month = factor(month(future_years)))
# Create seasonal dummy variables
future_data$month <- factor(month(future_data$index))
# Get the season
future_data$season <- case_when(
month(future_data$index) %in% c(12, 1, 2) ~ "Winter",
month(future_data$index) %in% c(3, 4, 5) ~ "Spring",
month(future_data$index) %in% c(6, 7, 8) ~ "Summer",
month(future_data$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Predict using the model
future_data$forecast <- predict(poly_month, newdata = future_data)
# Plot the forecasts
forecast_plot <- ggplot() +
geom_line(data = co2_tsib, aes(x = index, y = value), color = "black") +
geom_line(data = future_data, aes(x = index, y = forecast), color = "red") +
labs(title = "Forecast of Co2 Concentrations to 2020", x = "Year", y = "Co2 concentrations (ppmv)")
forecast_plot
# Generate future time points (e.g., monthly until 2020)
future_years <- seq(from = max(co2_tsib$index_date) + months(1), to = as.Date("2020-12-01"), by = "month")
future_data <- data.frame(index = future_years, month = factor(month(future_years)))
# Create seasonal dummy variables
future_data$month <- factor(month(future_data$index))
future_data$year <- factor(year(future_data$index))
# Get the season
future_data$season <- case_when(
month(future_data$index) %in% c(12, 1, 2) ~ "Winter",
month(future_data$index) %in% c(3, 4, 5) ~ "Spring",
month(future_data$index) %in% c(6, 7, 8) ~ "Summer",
month(future_data$index) %in% c(9, 10, 11) ~ "Autumn"
)
# Predict using the model
future_data$forecast <- predict(poly_month, newdata = future_data)
# Plot the forecasts
forecast_plot <- ggplot() +
geom_line(data = co2_tsib, aes(x = index, y = value), color = "black") +
geom_line(data = future_data, aes(x = index, y = forecast), color = "red") +
labs(title = "Forecast of Co2 Concentrations to 2020", x = "Year", y = "Co2 concentrations (ppmv)")
forecast_plot
# monthly time series with the line of best fit
co2_trend_plot <- co2_tsib %>%
ggplot(aes(x = index, y = value)) +
geom_line(color = 'black', size = .5) +
geom_smooth(method = "lm", formula = "y ~ x",se = F, color = 'blue', size = .8) +
labs(title = 'Average Trend Plotted on Monthly Average Co2 ppmv') +
xlab('Year') +
ylab('Co2 concentrations (ppmv)')
# average_yearly_increase
co2_tsib_yearly_change <- co2_tsib %>% as_tibble() %>%
mutate(year = year(index)) %>%
group_by(year) %>%
summarise(`yearly_co2` = mean(value)) %>%
ungroup() %>%
mutate(lag_co2 = lag(yearly_co2),
change = yearly_co2 - lag_co2,
percent_change = ((yearly_co2 - lag_co2)/yearly_co2)*100)
# getting average increase (i.e. size of the trend)
yearly_mean <- mean(co2_tsib_yearly_change$change, na.rm = T) # average 1.26 units of co2 change each year
yearly_sd <- sd(co2_tsib_yearly_change$change, na.rm = T) # with a sd of .51 units of co2
change_hist <- ggplot(co2_tsib_yearly_change, aes(x = change)) +
geom_histogram(color = 'gray20', fill = 'gray', binwidth = .25) +
scale_x_continuous(breaks = seq(floor(min(co2_tsib_yearly_change$change, na.rm =T)),
ceiling(max(co2_tsib_yearly_change$change, na.rm =T)), by = 0.5)) +
ggtitle('Histogram of Yearly changes in Co2 ppmv')
# bulk of increases yearly seem to be between .5 and 2 units on the
co2_trend_plot / change_hist