# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# +
# Imports
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from statsmodels.tsa.arima.model import ARIMA
sns.set_theme()
# Load data (city)
df = pd.read_csv('../data/GlobalLandTemperaturesByCity.csv')
df['dt'] = pd.DatetimeIndex(df['dt'])
df['Year'] = pd.DatetimeIndex(df['dt']).year
df = df[df['Year'] != 1743]
df = df[df['Year'] != 2013]

# For each country average all AverageTemperature values and add iso_alpha to countries
df_avg = pd.DataFrame({'AverageTemperature': df.groupby('Country')['AverageTemperature'].mean(), 'Country': df.groupby('Country')['Country'].first()})
df_avg = df_avg.reset_index(drop=True)
#print('Average country temperature df: ', df_avg)

df_avg_city = pd.DataFrame({'AverageTemperature': df.groupby('City')['AverageTemperature'].mean(), 'City': df.groupby('City')['City'].first()})
df_avg_city = df_avg_city.reset_index(drop=True)
#print('\n \n Average city temperature df: ', df_avg_city)

# For each Country for each Year average all AverageTemperature values
D1 = df.groupby(['Year', 'Country'])['AverageTemperature'].mean().reset_index()
meantemp = D1.groupby('Year')['AverageTemperature'].mean().reset_index() # Rough - er ren average.
meantemp1900 = meantemp[156:] # 1900-2012

# Normalize temperature relative to median
mediantemp = D1.groupby('Country')['AverageTemperature'].median()
df.groupby(['Year','Country'])['AverageTemperature'].mean().reset_index()['AverageTemperature']
#temp.columns = ['Year', 'Country', 'mean','median']
#D1['NormalizedTemperature'] = temp['mean'] - temp['median']
df

import copy
D2 = copy.deepcopy(D1)
normies = np.array([])
for country in D2['Country'].unique():
    norm = np.array(D2.groupby('Country')['AverageTemperature'].get_group(country) - D2.groupby('Country')['AverageTemperature'].get_group(country).median())
    normies = np.append(normies, norm)
D2['NormalizedTemperature'] = normies

# -

# Imports
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import plotly.express as px
import plotly.graph_objects as go
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_pacf, plot_acf
from statsmodels.tsa.stattools import adfuller
from itertools import product
from tqdm import tqdm_notebook
import warnings
warnings.filterwarnings("ignore")
sns.set_theme()

# +
# fit an ARIMA model and plot residual error

# fit model
modeltype = ARIMA(meantemp1900.AverageTemperature, order=(2,1,2))
model = modeltype.fit()
# summary of fit model
print(model.summary())
# line plot of residuals
residuals = pd.DataFrame(model.resid)
residuals.plot()

# # density plot of residuals
residuals.plot(kind='kde')
# # summary stats of residuals
print(residuals.describe())
# -

plt.plot(meantemp1900.AverageTemperature[2:])
plt.plot(model.fittedvalues[2:], color='red')


# ### Using TSA to predict future CO2 emissions

df3 = pd.read_csv('../data/annual-co-emissions-by-region.csv')
df3 = df3.rename(columns={'Annual CO2 emissions (zero filled)': "Annual CO2 emissions"})
df3 = df3.loc[df3['Entity'] == 'World']
df3 = df3.sort_values('Year')
px.line(df3, x = 'Year', y = 'Annual CO2 emissions')

# #### First order difference to account for increasing trend

diff1 = np.diff(np.log(df3['Annual CO2 emissions'])) 
df3['Log transformed first order differenced CO2 emissions'] = np.pad(diff1, (0, 1), 'constant')
df3['Log transformed CO2 emissions'] = np.log(df3['Annual CO2 emissions'])
px.line(df3, x = 'Year', y = 'Log transformed first order differenced CO2 emissions')

# #### Checking for stationarity

# Dickey-Fuller test of stationarity
ad_fuller_result = adfuller(df3['Log transformed first order differenced CO2 emissions'])
print(f'ADF Statistic: {ad_fuller_result[0]}')
print(f'p-value: {ad_fuller_result[1]}')

# #### Inspecting the ACF and PACF

acf = plot_acf(df3['Log transformed first order differenced CO2 emissions'], lags = 20)
pacf = plot_pacf(df3['Log transformed first order differenced CO2 emissions'], lags = 20)
plt.show(acf)
plt.show(pacf)


# Finding optimal ARIMA order
def optimize_ARIMA(order_list, exog):
    """
        Return dataframe with parameters and corresponding AIC
        
        order_list - list with (p, d, q) tuples
        exog - the exogenous variable
    """
    
    results = []
    
    for order in tqdm_notebook(order_list):
        try: 
            model = ARIMA(df3['Log transformed first order differenced CO2 emissions'], order = order).fit()
        except:
            continue
            
        aic = model.aic
        results.append([order, model.aic])
        
    result_df = pd.DataFrame(results)
    print(result_df)
    result_df.columns = ['(p, d, q)', 'AIC']
    #Sort in ascending order, lower AIC is better
    result_df = result_df.sort_values(by='AIC', ascending=True).reset_index(drop=True)
    
    return result_df


# +
ps = range(0, 8, 1)
d = 1
qs = range(0, 8, 1)# Create a list with all possible combination of parameters
parameters = product(ps, qs)
parameters_list = list(parameters)
order_list = []

for each in parameters_list:
    each = list(each)
    each.insert(1, 1)
    each = tuple(each)
    order_list.append(each)
    
result_df = optimize_ARIMA(order_list, exog=df3['Log transformed CO2 emissions'])
# -

print(result_df[result_df.AIC == result_df.AIC.min()])

# 1,1,2 ARIMA Model
model = ARIMA(df3['Log transformed CO2 emissions'], order=(1,1,2))
model_fit = model.fit()
print(model_fit.summary())
model_fit

y = 2020
for i in range(10):
    y = y + 1
    df3 = df3.append({'Entity': 'World', 'Code': 'OWID_WRL','Year': y, 'Annual CO2 emissions': float("NAN"), 'Log transformed first order differenced CO2 emissions': float("NAN"),'Log transformed CO2 emissions': float("NAN")}, ignore_index=True)

preds = model_fit.get_prediction(0,281)  # 95% conf
preds_ci = preds.conf_int()
preds_ci

preds_mu = preds.predicted_mean
df3['Predictions'] = preds_mu

# +
# Plot 
fig = px.line(df3, x='Year', y='Annual CO2 emissions')

# Only thing I figured is - I could do this 
fig.add_scatter(x=df3['Year'], y=np.exp(df3['Predictions']), mode='lines')

# Show plot 
fig.show()
# -



StarcoderdataJupyterScriptsDedupFiltered