From https://www.esrl.noaa.gov/gmd/ccgg/trends/data.html

Using Monthly Mean and Historical COmparison. Use Interpreted column for CO2 data.

Mauna Loa CO2 data remove text leaving 1 header row. Ranamed colums.

library(forecast)

Read

mlmco2 <- read.csv(file = 'co2_mm_mlo.csv')
cat('Min Date:',min(mlmco2$year),' Max date:',max(mlmco2$year))
## Min Date: 1958  Max date: 2019

Exploratory Plots

Note Strong Tailing at the low and high end.

mlm21<-mlmco2['interp']

par( mfrow=c(1,2) )
hist(mlm21$interp, main='Mauna Loa CO2 Monthly', xlab='CO2 PPM')
qqnorm(mlm21$interp,main='Mauna Loa CO2 Monthly \n NORM Q-Q Plot')
qqline(mlm21$interp)

Create Time Series

Create time series and make plots

ts.mlm21 <-ts(mlm21$interp,frequency = 12, start=c(1958,3))
par( mfrow=c(1,1) )
plot.ts(ts.mlm21, main='Mauna Loa CO2 (PPM)', ylab='PPM')

Note slow decay characteristic of seasonal processes.

acf(ts.mlm21,main='Mauna Loa CO2 PPM - AR')

Note seasonal repition of significant lags.

pacf(ts.mlm21,main='Mauna Loa CO2 PPM - MA')

Auto ARMIA Analysis

Indicates seasonality observed

ts.mlm21.asa <- auto.arima(ts.mlm21)
ts.mlm21.asa
## Series: ts.mlm21 
## ARIMA(1,1,1)(2,1,2)[12] 
## 
## Coefficients:
##          ar1      ma1     sar1     sar2     sma1     sma2
##       0.1920  -0.5626  -0.3679  -0.0303  -0.4944  -0.3120
## s.e.  0.0903   0.0765   0.7052   0.0435   0.7048   0.6099
## 
## sigma^2 estimated as 0.09853:  log likelihood=-188.41
## AIC=390.83   AICc=390.98   BIC=422.97

Holt Winters to Evaluate Seasonality

hw.ts.mlm21<-HoltWinters(ts.mlm21)
hw.ts.mlm21
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = ts.mlm21)
## 
## Smoothing parameters:
##  alpha: 0.5157221
##  beta : 0.01470134
##  gamma: 0.3155293
## 
## Coefficients:
##            [,1]
## a   412.1896475
## b     0.1958788
## s1    0.5909013
## s2    1.0600260
## s3    1.6925863
## s4    2.9995738
## s5    3.5685586
## s6    2.5196995
## s7    0.4517418
## s8   -1.6562657
## s9   -3.2470785
## s10  -3.1066365
## s11  -1.6322918
## s12  -0.4720103
hw.ts.mlm21$SSE
## [1] 74.56695

Obtain Forecast

fc.hw.ts.mlm21<-forecast(hw.ts.mlm21,25)
fc.hw.ts.mlm21
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2020       412.9764 412.5678 413.3851 412.3515 413.6014
## Feb 2020       413.6414 413.1802 414.1026 412.9361 414.3468
## Mar 2020       414.4699 413.9602 414.9795 413.6904 415.2493
## Apr 2020       415.9727 415.4176 416.5279 415.1237 416.8218
## May 2020       416.7376 416.1393 417.3359 415.8226 417.6526
## Jun 2020       415.8846 415.2450 416.5243 414.9064 416.8629
## Jul 2020       414.0125 413.3330 414.6921 412.9733 415.0518
## Aug 2020       412.1004 411.3822 412.8186 411.0021 413.1988
## Sep 2020       410.7055 409.9497 411.4613 409.5496 411.8614
## Oct 2020       411.0418 410.2492 411.8344 409.8297 412.2539
## Nov 2020       412.7120 411.8834 413.5406 411.4448 413.9793
## Dec 2020       414.0682 413.2042 414.9322 412.7468 415.3896
## Jan 2021       415.3270 414.4089 416.2451 413.9229 416.7311
## Feb 2021       415.9920 415.0402 416.9438 414.5364 417.4476
## Mar 2021       416.8204 415.8353 417.8055 415.3138 418.3270
## Apr 2021       418.3233 417.3051 419.3414 416.7662 419.8804
## May 2021       419.0881 418.0373 420.1390 417.4809 420.6953
## Jun 2021       418.2352 417.1518 419.3186 416.5782 419.8921
## Jul 2021       416.3631 415.2474 417.4788 414.6568 418.0694
## Aug 2021       414.4510 413.3031 415.5988 412.6955 416.2064
## Sep 2021       413.0560 411.8763 414.2358 411.2517 414.8603
## Oct 2021       413.3923 412.1808 414.6039 411.5394 415.2453
## Nov 2021       415.0626 413.8193 416.3058 413.1612 416.9640
## Dec 2021       416.4187 415.1439 417.6936 414.4690 418.3684
## Jan 2022       417.6775 416.3562 418.9989 415.6567 419.6984

Plot with Forecast

s <- 'Mauna Loa Co2 (PPM) Meas + ' 
sp <- 'Forecast'
par( mfrow=c(1,1) )
plot(mlmco2$decdata, mlmco2$interp, type='o',  pch=20, 
     main='Mauna Loa Co2 (PPM) Meas + Forecast',
     xlab='Year', ylab='PPM',
     col=ifelse( mlmco2$interp>=400,'blueviolet','blue'),
     xlim=c(1958,2022), ylim=c(310,420))
lines(fc.hw.ts.mlm21$mean, type='o', pch=20, col='red')
legend(1960, 420, legend=c("Meas < 400 PPM", "Mea >= 400 PPM", "Forecast"),
       fill=c("blue", "blueviolet", 'red'), cex=0.8, bg='lightblue')

END