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

# ### Predict CO2 emissions from Mauna Lao volcano

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.model_selection import train_test_split
from statsmodels.tsa.deterministic import DeterministicProcess
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression


co2_df = pd.read_csv("../ts-course-data/co2.csv", index_col = "Date", parse_dates=True)
co2_df.index.freq = "MS"

co2_df.head(5)

co2_df.plot(figsize = (15, 7), title = "Mauna Lao CO2 emissions")
plt.show()

# +
# It does look like there is a linear trend plus seasonality

# +
# fit a basic linear regression model using DeterministicProcess
# split the df into train and test, lets use 10% data as the test dataset

dp = DeterministicProcess (
    index=co2_df.index,  # dates from the training data
    constant=True,       # dummy feature for the bias (y_intercept)
    order=1,             # the time dummy (trend)
    drop=True
)

X = dp.in_sample()
y = co2_df.loc[:, "CO2"]

# -

X_train, X_test = train_test_split(X, test_size = 0.1, shuffle = False)
y_train, y_test = train_test_split(y, test_size = 0.1, shuffle = False)

X_train.head(4), y_train.head(4)

X_test.head(4), y_test.head(4)

# +
model = LinearRegression()
model.fit(X_train, y_train)

y_pred = pd.Series(model.predict(X_test), index=X_test.index)
# -

rmse = np.sqrt(mean_squared_error(y_pred, y_test))
print("rmse is {}".format(rmse))

# ### RMSE using basic linear regression model on just time dummy trend is 3.34

y_pred.head(4)

y_test.head(4)

plt.figure(figsize = (10, 5))
ax = plt.subplot(1, 1, 1)
ax.plot(y_test, label = "y test")
ax.plot(y_pred, color = "red", label = "y pred")
plt.legend()
plt.show()

# +
# Lets try to add seasonality.
model = LinearRegression()
model.fit(X, y)

y_pred = pd.Series(model.predict(X), index=X.index)
plt.figure(figsize = (10, 5))
ax = plt.subplot(1, 1, 1)
ax.plot(y - y_pred, color = "red", label = "Yactual - Ypred")
plt.title("Seasonality = Yactual - Ypred(includes trends)")
plt.legend()
plt.show()

# +
# Merging the trend + seasonality
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess

fourier = CalendarFourier(freq="A", order=10)  # 10 sin/cos pairs for "A"nnual seasonality

dp = DeterministicProcess(
    index=co2_df.index,
    constant=True,               # dummy feature for bias (y-intercept)
    order=1,                     # trend (order 1 means linear)
    seasonal=True,               # monthly seasonality (indicators)
    additional_terms=[fourier],  # annual seasonality (fourier)
    drop=True,                   # drop terms to avoid collinearity
)

X = dp.in_sample() 
y = co2_df.loc[:, "CO2"]

X.head(5)
# -

X_train, X_test = train_test_split(X, test_size = 0.1, shuffle = False)
y_train, y_test = train_test_split(y, test_size = 0.1, shuffle = False)

# +
model = LinearRegression()
model.fit(X_train, y_train)

y_pred = pd.Series(model.predict(X_test), index=X_test.index)
# -

rmse = np.sqrt(mean_squared_error(y_pred, y_test))
print("rmse is {}".format(rmse))

# ### RMSE is reduced from 3.34 (just linear trend) to 2.51 using seasonality (monthly and yearly)

y_pred.head(4)

y_test.head(5)

plt.figure(figsize = (10, 5))
ax = plt.subplot(1, 1, 1)
ax.plot(y_test, label = "y test")
ax.plot(y_pred, color = "red", label = "y pred")
plt.title("Y Test vs Y Pred")
plt.legend()
plt.show()

# +
# Lets try to add cycles.
model = LinearRegression()
model.fit(X, y)

y_pred = pd.Series(model.predict(X), index=X.index)
plt.figure(figsize = (10, 5))
ax = plt.subplot(1, 1, 1)
ax.plot(y - y_pred, color = "red", label = "Yactual - Ypred")
plt.title("Cycles = Yactual - Ypred(trends + seasonality)")
plt.legend()
plt.show()

# +
# lets add a lag of 1

from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess

fourier = CalendarFourier(freq="A", order=10)  # 10 sin/cos pairs for "A"nnual seasonality

dp = DeterministicProcess(
    index=co2_df.index,
    constant=True,               # dummy feature for bias (y-intercept)
    order=1,                     # trend (order 1 means linear)
    seasonal=True,               # weekly seasonality (indicators)
    additional_terms=[fourier],  # annual seasonality (fourier)
    drop=True,                   # drop terms to avoid collinearity
)

X = dp.in_sample() 
y = co2_df.loc[:, "CO2"]

X["lag_1"] = y.shift(1) # adding a lag of 1 for cycles
X = X.iloc[1: , :] # dropping the first row with NaN
y = y[1:]
X_train, X_test = train_test_split(X, test_size = 0.1, shuffle = False)
y_train, y_test = train_test_split(y, test_size = 0.1, shuffle = False)

model = LinearRegression()
model.fit(X_train, y_train)

y_pred = pd.Series(model.predict(X_test), index=X_test.index)

rmse = np.sqrt(mean_squared_error(y_pred, y_test))
print("rmse is {}".format(rmse))
# -

plt.figure(figsize = (10, 5))
ax = plt.subplot(1, 1, 1)
ax.plot(y_test, label = "y test")
ax.plot(y_pred, color = "red", label = "y pred")
plt.title("Y Test vs Y Pred")
plt.legend()
plt.show()

# +
# get the residual

model = LinearRegression()
model.fit(X, y)

y_pred = pd.Series(model.predict(X), index=X.index)
plt.figure(figsize = (10, 5))
ax = plt.subplot(1, 1, 1)
ax.plot(y - y_pred, color = "red", label = "Yactual - Ypred")
plt.title("Residual = Yactual - Ypred(trends + seasonality + cycles)")
plt.legend()
plt.show()

# +
# Train XGBoost on the residuals
from sklearn.ensemble import GradientBoostingRegressor
y_residual = y - y_pred

X_train_residual, X_test_residual = train_test_split(X, test_size = 0.1, shuffle = False)
y_train_residual, y_test_residual = train_test_split(y_residual, test_size = 0.1, shuffle = False)

model = GradientBoostingRegressor()
model.fit(X_train_residual, y_train_residual)

y_pred_residual = pd.Series(model.predict(X_test_residual), index=X_test_residual.index)

rmse = np.sqrt(mean_squared_error(y_pred_residual, y_test_residual))
print("rmse is {}".format(rmse))
# -

y_test_residual.head(4), y_pred_residual.head(4)

len(y_test_residual), len(y_pred_residual)

# +
model = GradientBoostingRegressor()
model.fit(X, y_residual)

y_pred_residual = pd.Series(model.predict(X), index=X.index)
# -

# overall y is the sum of y predicted from trends + seasonality + cycles (from Linear Regression) and 
# y predicted from residual using XGBoostRegressor
final_y = y_pred + y_pred_residual

print(y_pred.head(4), y_pred_residual.head(4), final_y.head(4))

y_pred, final_y

plt.figure(figsize = (15, 6))
ax = plt.subplot(1, 1, 1)
ax.plot(final_y, color = "red", label = "Y Pred", alpha = 0.5)
ax.plot(y, color = "blue", label = "Actual Y", alpha = 0.2)
plt.title("Final model - Linear + GradientBoost")
plt.show()

# +
_, final_y_test = train_test_split(final_y, test_size = 0.1, shuffle = False)

rmse = np.sqrt(mean_squared_error(final_y_test, y_test))
print("rmse is {}".format(rmse))
# -

# ### RMSE using combination of both Linear and Gradient Boost is 0.149 on the test dataset

plt.figure(figsize = (10, 5))
ax = plt.subplot(1, 1, 1)
ax.plot(y_test, label = "y test")
ax.plot(final_y_test, color = "red", label = "y pred")
plt.title("Y Test vs Y Pred")
plt.legend()
plt.show()


# +
# Perfecto!

StarcoderdataJupyterScriptsDedupFiltered