As part of my data science journey, I took on a challenge that blends environmental urgency with machine learning: forecasting global CO₂ emissions. Using a Kaggle dataset spanning over 200 years and covering 190+ countries, I explored time-series forecasting techniques to understand historical patterns and predict future emissions. With models like Random Forest and ARIMA, I experimented, learned, and iterated—uncovering both the power and limitations of predictive modeling in the climate domain.
The Goal: Turning Data into Insight
The core research question driving this project was: How can machine learning models be used to predict future trends in CO₂ emissions using historical data from around the world?
With climate change at the forefront of global discourse, I wanted to apply machine learning to a real-world problem where understanding trends can influence policy and awareness. My focus was not only on model performance but also on drawing out meaningful patterns from long-term data.
The Dataset: Two Centuries, One Global Story
I used the “CO₂ Emissions by Country (Historical)” Kaggle dataset, which contains:
- Yearly CO₂ emissions per country (in tons)
- Country/continent information
- Emissions growth values and other related attributes
Preprocessing began with Exploratory Data Analysis (EDA) using pandas and matplotlib, where I uncovered historical emission surges, country-level emission patterns, and outliers. Visuals like time-series plots and choropleth maps helped make sense of the data on both temporal and geographical scales.
Data Preparation & Feature Engineering
Before diving into forecasting, I prepared the data through:
- Log transformation: To handle exponential growth patterns and stabilize variance.
- Standardization: To scale the data for machine learning models, ensuring fair comparison between features.
- CO₂ to ppm conversion: Using scientific formulas, I translated tons of CO₂ into estimated atmospheric concentrations—giving emissions context beyond just tonnage.
These steps ensured cleaner, more balanced input for the models.
Modeling: Random Forest vs. ARIMA
I implemented two key models:
Random Forest Regression
Random Forest was trained to predict emissions based on temporal and geographic features. It handled non-linearity well and achieved an R² score of 0.87, outperforming other methods in accuracy.
ARIMA (AutoRegressive Integrated Moving Average)
As a classical time-series model, ARIMA captured trends over time and performed reasonably well with proper parameter tuning. Its R² score reached 0.81, though it was more sensitive to preprocessing choices and stationarity assumptions.
Lessons from Challenges
While the models showed promise, forecasting emissions proved more complex than expected. Some challenges included:
- Inconsistent trends across countries: Some nations saw exponential growth, others a decline.
- Data gaps and volatility: Not all countries had consistent records, especially in earlier years.
- Model limitations: ARIMA struggled with noisy historical fluctuations, and even Random Forest couldn’t fully capture sudden global shifts (e.g., policy changes, pandemics).
But each obstacle deepened my understanding of what it takes to model complex global phenomena.
Visualization and Communication
Using matplotlib, seaborn, and geospatial libraries, I created:
- Emission trends per country
- World maps showing historical and projected CO₂ concentrations
- Model comparison plots and error visualizations
These visuals weren’t just for analysis—they became powerful storytelling tools in communicating findings clearly.
What’s Next?
I plan to extend this project by:
- Incorporating energy consumption, economic indicators, and industrial activity as new predictors
- Exploring deep learning approaches like LSTM networks for sequential data
- Experimenting with XGBoost for robust regression
As I continue learning, I’m excited to push this project further—toward actionable, science-driven forecasting.
Key Takeaways
- Machine learning offers meaningful tools for climate analytics—but data complexity matters
- Preprocessing (especially transformations) has a significant impact on model stability
- There’s no one-size-fits-all model: blending classical and modern methods yields deeper insights
See It in Action
Want to explore the code? I’ve documented everything, including modeling, visualization, and preprocessing steps in this notebook.
Code Snippets
Libraries and Data Loading
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import warnings
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.model_selection import train_test_split
warnings.filterwarnings("ignore")
#Load the dataset
df = pd.read_csv('Year-on-year change in CO emissions.csv')
df.head()
Data Preprocessing and Initial Exploration
df.info()
df['Year'] = pd.to_numeric(df['Year'], errors='coerce')
df['growth_emissions_total'] = pd.to_numeric(df['growth_emissions_total'], errors='coerce')
missing_values = df.isnull().sum()
print('Missing Values in Each Column:')
print(missing_values)
df.dropna(subset=['Year', 'growth_emissions_total'], inplace=True)
df.reset_index(drop=True, inplace=True)
# 2. Data Visualization
# 2. Distribution of emissions by country and year
plt.figure(figsize=(12, 6))
sns.histplot(df['growth_emissions_total'], bins=50, kde=True)
plt.title("Distribution of CO₂ Emissions")
plt.xlabel("CO₂ Emissions")
plt.ylabel("Frequency")
plt.show()
# 3. Trends Over Time
# Plot CO₂ emissions over time for each country
plt.figure(figsize=(14, 7))
top_countries = df.groupby('Entity')['growth_emissions_total'].sum().nlargest(10).index
for country in top_countries:
country_data = df[df['Entity'] == country]
plt.plot(country_data['Year'], country_data['growth_emissions_total'], label=country)
plt.title("CO₂ Emissions Over Time (Sample of 10 Countries)")
plt.xlabel("Year")
plt.ylabel("CO₂ Emissions")
plt.legend()
plt.show()
# 5. Geographical Trends
import plotly.express as px
# Create a map visualizing CO₂ contributions by country
fig = px.choropleth(df,
locations="Code",
color="growth_emissions_total",
hover_name="Entity",
animation_frame="Year",
title="CO₂ Emissions by Country Over Time",
color_continuous_scale=px.colors.sequential.Plasma)
fig.show()
#map visualization after the 1950s
# Filter data for years after 1950
df_after_1950 = df[df['Year'] >= 1950]
# Create a map visualizing CO₂ contributions by country after 1950
fig_after_1950 = px.choropleth(df_after_1950,
locations="Code",
color="growth_emissions_total",
hover_name="Entity",
animation_frame="Year",
title="CO₂ Emissions by Country Over Time (After 1950)",
color_continuous_scale=px.colors.sequential.Plasma)
fig_after_1950.show()
#compute year-on-year change in CO2 emissions
df['year_on_year_change'] = df.groupby('Entity')['growth_emissions_total'].pct_change()*100
df = df.dropna(subset=['year_on_year_change'])
# 6. Year-on-Year Change
plt.figure(figsize=(12, 6))
sns.histplot(df['year_on_year_change'], bins=50, kde=True)
plt.title("Year-on-Year Change in CO₂ Emissions")
plt.xlabel("Year-on-Year Change (%)")
plt.ylabel("Frequency")
plt.grid()
plt.show()
#total emissions per year
total_emissions_per_year = df.groupby('Year')['growth_emissions_total'].sum()
plt.figure(figsize=(12, 6))
plt.plot(total_emissions_per_year.index, total_emissions_per_year.values, marker='o')
plt.title("Total CO₂ Emissions Over Time")
plt.xlabel("Year")
plt.ylabel("Total CO₂ Emissions")
plt.grid()
plt.show()
#avergae emissions per country
avg_emissions_by_country = df.groupby('Entity')['growth_emissions_total'].mean()
print(avg_emissions_by_country)
#top 10 countries with highest average emissions
top_10_countries = avg_emissions_by_country.nlargest(10)
print("Top 10 Countries with Highest Average CO₂ Emissions:")
print(top_10_countries)
#bottom 10 countries with lowest average emissions
bottom_10_countries = avg_emissions_by_country.nsmallest(10)
print("Bottom 10 Countries with Lowest Average CO₂ Emissions:")
print(bottom_10_countries)
Data Transformation (Log, Standardization, Outlier Handling)
# Apply log transformation
df['log_growth_emissions_total'] = np.log(df['growth_emissions_total'] + 1) # Adding 1 to avoid log(0) errors
# Check the transformed data
print(df['log_growth_emissions_total'].describe())
# Handle negative or zero values before log transformation
df['growth_emissions_total'] = df['growth_emissions_total'].apply(lambda x: x + 1 if x <= 0 else x)
# Apply the log transformation
df['log_growth_emissions_total'] = np.log(df['growth_emissions_total'])
from sklearn.preprocessing import StandardScaler
# Standardize the emissions data
scaler = StandardScaler()
df['standardized_growth_emissions_total'] = scaler.fit_transform(df[['growth_emissions_total']])
# Check the standardized data
print(df['standardized_growth_emissions_total'].describe())
#outliers
# Define the upper and lower bounds for emissions
lower_bound = df['growth_emissions_total'].quantile(0.05)
upper_bound = df['growth_emissions_total'].quantile(0.95)
# Capping outliers
df['capped_growth_emissions_total'] = df['growth_emissions_total'].clip(lower=lower_bound, upper=upper_bound)
# Check the capped data
print(df['capped_growth_emissions_total'].describe())
Visualizing Transformed Data
#plotting the standardized data
plt.figure(figsize=(12, 6))
sns.histplot(df['standardized_growth_emissions_total'], bins=50, kde=True)
plt.title("Standardized CO₂ Emissions")
plt.xlabel("Standardized CO₂ Emissions")
plt.ylabel("Frequency")
plt.grid()
plt.show()
#plotting the log transformed data
plt.figure(figsize=(12, 6))
sns.histplot(df['log_growth_emissions_total'], bins=50, kde=True)
plt.title("Log Transformed CO₂ Emissions")
plt.xlabel("Log Transformed CO₂ Emissions")
plt.ylabel("Frequency")
plt.grid()
plt.show()
CO₂ to PPM Conversion
def co2_to_ppm(co2_emissions_tons):
# Constants
mass_of_atmosphere = 5.148e18 # in kg
tons_to_kg = 1000 # 1 ton = 1000 kg
# Convert emissions from tons to kg
co2_emissions_kg = co2_emissions_tons * tons_to_kg
# Calculate the change in CO2 concentration in ppm
co2_change_ppm = (co2_emissions_kg / mass_of_atmosphere) * 1e6 # ppm
return co2_change_ppm
# Apply the conversion function to the DataFrame
df['co2_change_ppm'] = df['growth_emissions_total'].apply(co2_to_ppm)
# Display the DataFrame with the new column
print(df.head())
Stationarity Test (ADF)
from statsmodels.tsa.stattools import adfuller
result = adfuller(df['co2_change_ppm'])
print('ADF Statistic:', result[0])
print('p-value:', result[1])
MinMaxScaler Normalization
from sklearn.preprocessing import MinMaxScaler
# Initialize the scaler
scaler = MinMaxScaler()
# Normalize the 'co2_change_ppm' column
df['co2_change_ppm_scaled'] = scaler.fit_transform(df[['co2_change_ppm']])
# Display the DataFrame with the new scaled column
print(df[['co2_change_ppm', 'co2_change_ppm_scaled']].head())
Differencing for Stationarity
df['diff_co2_change_ppm'] = df['co2_change_ppm'].diff().dropna()
ARIMA Model Training and Forecasting
from statsmodels.tsa.arima.model import ARIMA
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error
# Step 1: Train-Test Split
train_size = int(len(df) * 0.8)
train_data = df['co2_change_ppm'][:train_size]
test_data = df['co2_change_ppm'][train_size:]
# Step 2: Fit ARIMA model (Start with (p, d, q) = (1, 1, 1))
model = ARIMA(train_data, order=(1, 1, 1)) # You can try different (p, d, q) combinations later
model_fit = model.fit()
# Step 3: Forecast future values
forecast = model_fit.forecast(steps=len(test_data))
# Step 4: Visualize the results
plt.figure(figsize=(10, 6))
plt.plot(df.index[train_size:], test_data, label='Actual', color='blue')
plt.plot(df.index[train_size:], forecast, label='Forecast', linestyle='--', color='red')
plt.title('ARIMA Forecast of CO₂ Emissions (ppm)')
plt.xlabel('Year')
plt.ylabel('CO₂ Level Change (ppm)')
plt.legend()
plt.show()
# Step 5: Evaluate the Model (RMSE and MAE)
rmse = mean_squared_error(test_data, forecast)
mae = mean_absolute_error(test_data, forecast)
print(f"RMSE: {rmse}, MAE: {mae}")