
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
%matplotlib inline

from statsmodels.tsa.seasonal import seasonal_decompose
from datetime import datetime, timedelta
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from random import gauss
from random import seed
import scipy.stats as stats
import math
import os
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
#from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.arima_model import ARIMA


# Style Dictionary to standarize plotting 
# scheme between different python scripts 
styledict = {'xtick.labelsize':21,
             'xtick.major.size':9,
             'xtick.major.width':1,
             'ytick.labelsize':21,
             'ytick.major.size':9,
             'ytick.major.width':1,
             'legend.framealpha':0.0,
             'legend.fontsize':15,
             'axes.labelsize':22,
             'axes.titlesize':22,
             'axes.linewidth':2,
             'figure.figsize':(12,8),
             'savefig.format':'jpg'}
plt.rcParams.update(**styledict)

# define font style for text
font = {'family': 'serif',
        'color':  'blue',
        'weight': 'normal',
        'size': 12,
        }


# Go to working directory 
os.chdir("/users/sshen/climstats")


#Python code for plotting: Fig. 7.1: The Keeling curve
os.chdir("/users/sshen/climstats")
co2m = pd.read_csv('data/co2_mm_mlo.csv', skiprows = 51)
print(co2m.head(5))
df = pd.DataFrame(co2m)
print(df.shape)
plt.plot(co2m['decimal date'], co2m['average'], color = 'r')
plt.plot(co2m['decimal date'], co2m['interpolated'], 
         color = 'k')
plt.title('Monthly Atmospheric CO\N{SUBSCRIPT TWO} \
at Mauna Loa Observatory', pad = 15)
plt.ylabel('Parts Per Million [ppm]',
          labelpad = 15)
# text inside the figure
plt.text(1958,410, 
      s='Scripps Institution of Oceanography', size = 18)
plt.text(1958,400, 
         s='NOAA Global Monitoring Laboratory', size = 18)
# save figure
plt.savefig("fig0701.eps")

   year  month  decimal date  average  interpolated  trend  ndays  Unnamed: 7
0  1958      3     1958.2027   315.70        314.43     -1  -9.99       -0.99
1  1958      4     1958.2877   317.45        315.16     -1  -9.99       -0.99
2  1958      5     1958.3699   317.51        314.71     -1  -9.99       -0.99
3  1958      6     1958.4548   317.24        315.14     -1  -9.99       -0.99
4  1958      7     1958.5370   315.86        315.18     -1  -9.99       -0.99
(749, 8)


os.chdir("/users/sshen/climstats")
co2m = pd.read_csv('data/co2_mm_mlo.csv', skiprows = 51)


# get date into correct format in order to plot time series: 
date_index = []

# convert decimal date into datetime object, and store
df_co2 = df 
for i in df_co2['decimal date']:
    start = float(i)
    year = int(float(start))
    rem = start - year
    base = datetime(year, 1, 1)
    result = base + timedelta(seconds=(base.replace(year=
           base.year + 1) - base).total_seconds() * rem)
    date_index.append(result)

date_index = pd.to_datetime(date_index)


# Create a time series starting from March 1958 to July 2020 with monthly frequency
from statsmodels.tsa.seasonal import seasonal_decompose
os.chdir("/users/sshen/climstats")
co2m = pd.read_csv('data/co2_mm_mlo.csv', skiprows = 51) 
data_values = ser = pd.Series(co2m['average'])
ser.index = date_index
# get observed data in correct format
# convert values from string to float):
data_values = [float(i) for i in ser]
co2 = ser
co2.index = pd.date_range('1958-03-17', '2020-08-17', freq = 'M')
#pd.date_range(start='1958-03', end='2020-08', freq='M')

# Display the time series
print(co2.head())

# Decompose the time series into trend, seasonal, and residual components
decomposition = sm.tsa.seasonal_decompose(co2, model='additive')

# Plot the decomposed time series
plt.figure(figsize=(11, 16))
plt.subplot(411)
plt.plot(co2)
plt.title(r'Monthly Atmospheric CO$_2$ [ppm] at Mauna Loa Observatory')
plt.ylabel('Observed Data')
plt.subplot(412)
plt.plot(decomposition.trend)
plt.ylabel('Trend')
plt.subplot(413)
plt.plot(decomposition.seasonal)
plt.ylabel('Seasonality')
plt.subplot(414)
plt.plot(decomposition.resid)
plt.ylabel('Random')
plt.tight_layout()
plt.savefig("fig0702.eps")
plt.show()

#you can also use a simple command to make the plots
#decomposition.plot()
#but the figures do not look as nice

1958-03-31    315.70
1958-04-30    317.45
1958-05-31    317.51
1958-06-30    317.24
1958-07-31    315.86
Freq: M, Name: average, dtype: float64

Fig. 7.1: Monthly CO2 Keeling Curve¶

Fig. 7.2: Monthly Atmospheric CO2 at Mauna Loa Observatory¶