Luis José Zapata Blog

Python basics

Fri, 01 Jan 2021 00:00:00 +0000

Build a foundation in Python.

1-2 hours per week, for 8 weeks

Learn

Quiz

What is the difference between lists and tuples?

Lists

Lists are mutable - they can be changed
Slower than tuples
Syntax: a_list = [1, 2.0, 'Hello world']

Tuples

Tuples are immutable - they can’t be changed
Tuples are faster than lists
Syntax: a_tuple = (1, 2.0, 'Hello world')

Is Python case-sensitive?

Yes

Visualization

Fri, 01 Jan 2021 00:00:00 +0000

Learn how to visualize data with Plotly.

1-2 hours per week, for 8 weeks

Learn

Quiz

When is a heatmap useful?

Lorem ipsum dolor sit amet, consectetur adipiscing elit.

Write Plotly code to render a bar chart

import plotly.express as px
data_canada = px.data.gapminder().query("country == 'Canada'")
fig = px.bar(data_canada, x='year', y='pop')
fig.show()

Statistics

Fri, 01 Jan 2021 00:00:00 +0000

Introduction to statistics for data science.

1-2 hours per week, for 8 weeks

Learn

The general form of the normal probability density function is:

$$ f(x) = \frac{1}{\sigma \sqrt{2\pi} } e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} $$

The parameter $\mu$ is the mean or expectation of the distribution. $\sigma$ is its standard deviation. The variance of the distribution is $\sigma^{2}$.

Quiz

What is the parameter $\mu$?

The parameter $\mu$ is the mean or expectation of the distribution.

Using Python, SerpAPI, and GPT to Build a Smart Web Search Analyzer

Wed, 12 Mar 2025 00:00:00 +0000

Introduction

In this blog post, I’ll walk through how to use modern tools like Python, SerpAPI, and OpenAI’s GPT models to build a modular search pipeline. The idea is to scrape Google Search results programmatically, analyze them using LLMs, and extract structured insights.

What Are APIs, and Why Are We Using Them?

APIs (Application Programming Interfaces) are like bridges that allow one software to talk to another in a standardized way. Instead of building a search engine or a language model from scratch, developers can connect to existing services using APIs to request data, process information, or trigger tasks.

In this project, we rely on two APIs:

SerpAPI: This is a powerful tool that lets us programmatically fetch Google search results. Instead of scraping raw HTML and fighting with unpredictable webpage structures, SerpAPI gives us clean, structured JSON data for any Google query — including titles, URLs, and summaries.
OpenAI API: This allows us to access large language models (like GPT-4o) via code. By sending a structured prompt and some example text to OpenAI’s API, we receive intelligent summaries, classifications, or decisions in return — perfect for analyzing search results at scale.

Together, these APIs let us build a small but powerful pipeline: search → filter → analyze — all from Python.

To keep things simple and general, our example will revolve around analyzing search results related to shoes, rather than any specialized domain, but this project can easily be adapted to any other complex topic.

We’ll break the process into four key parts:

Fetching Google results using SerpAPI
Cleaning and organizing the raw data
Using GPT to extract structured insight from results
Putting everything together in one function

Step 1: Fetch Google Results Using SerpAPI

We’ll start by querying SerpAPI. SerpAPI is a commercial API that allows developers to access real-time, structured Google search results without the hassle of traditional web scraping. In this blog, SerpAPI is used to automate the collection of search result data for a given keyword or product — for example, “shoes with waterproof soles” or “trail running sneakers.” Instead of manually performing searches and copying links, we send queries through SerpAPI and receive a JSON response containing the search titles, URLs, summaries, and even metadata like the publication date. This structured output is the foundation of our analysis pipeline, feeding directly into the next steps where we clean the data and evaluate relevance using language models.

First, we’ll need a SerpAPI API key.

import requests
def fetch_google_results(query, api_key, start=0, num=10):
params = {
"q": query,
"location": "Los Angeles, California, United States",
"start": start,
"num": num,
"hl": "en",
"gl": "us",
"engine": "google",
"api_key": api_key
}
response = requests.get("https://serpapi.com/search", params=params)
results = response.json().get("organic_results", [])
return results

Example:

import pandas as pd
serpapi_key = serapi_code
results = fetch_google_results("best running shoes 2024", serpapi_key)
print(pd.DataFrame(results).head())

## position ... rich_snippet
## 0 1 ... NaN
## 1 2 ... NaN
## 2 3 ... NaN
## 3 4 ... NaN
## 4 5 ... NaN
##
## [5 rows x 13 columns]

Keep in mind that SerpAPI offers 100 free requests per month. If your project requires more extensive usage, you’ll need to upgrade to a paid plan. Additionally, SerpAPI supports other search engines (like Bing, DuckDuckGo, YouTube, and more), so you can adapt the same API logic to different platforms if needed.

Step 2: Clean and Organize the Raw Data

Once we have the raw search results, we want to extract important fields and handle possible issues like missing data.

import pandas as pd
import datetime
def clean_results(raw_results):
df = pd.DataFrame(raw_results)
df = df[["position", "title", "link", "snippet"]].rename(
columns={
"position": "Rank",
"title": "Title",
"link": "URL",
"snippet": "Summary"
}
)
df["Scrape Date"] = datetime.date.today().strftime("%Y-%m-%d")
return df

Example:

df_cleaned = clean_results(results)
df_cleaned.head()

## Rank ... Scrape Date
## 0 1 ... 2025-04-12
## 1 2 ... 2025-04-12
## 2 3 ... 2025-04-12
## 3 4 ... 2025-04-12
## 4 5 ... 2025-04-12
##
## [5 rows x 5 columns]

Step 3: Use GPT to Analyze the Results

OpenAI’s API provides access to advanced language models like GPT-4, which can process and understand natural language. In this project, we use the Chat Completion API to enhance our search pipeline by analyzing each result retrieved from Google. For every link found, we prompt GPT to generate a brief summary of the content and assign a probability score that estimates how likely the link is to contain useful shopping-related information — such as detailed product reviews, buying guides, or comparisons. This allows us to automatically prioritize the most informative links without manual filtering, making it easier to identify high-value content in large batches of search results.

from openai import OpenAI
import json
import time
client = OpenAI(api_key=openai_code)
def query_openai(df, keyword):
json_data = df[["Rank", "Title", "URL", "Summary"]].to_dict(orient="records")
prompt = f""" You are an AI assistant helping to analyze search engine results for consumer products. Analyze the following JSON search results from Google. For each result, please:
1. Provide a short summary (max 150 characters) of the result.
2. Estimate the probability (from 0 to 1) that the result contains useful information such as reviews, comparisons, product specs, or guides relevant to the keyword '{keyword}'.
Return a JSON list under the key "results" with these fields:
- Rank
- URL
- Probability
- GPT_Summary
Ensure that the number of output rows matches the number of input rows. If any row has no relevant information, still return it with a low Probability.
####
Data:
{json.dumps(json_data, indent=2)}
####
"""
response = client.chat.completions.create(
model="gpt-4o",
messages=[
{"role": "system", "content": "You are an AI summarizer and classifier of search engine results."},
{"role": "user", "content": prompt}
],
response_format={
"type": "json_schema",
"json_schema": {
"name": "SearchEvaluation",
"schema": {
"type": "object",
"properties": {
"Rank": {"type": "integer"},
"URL": {"type": "string"},
"Probability": {"type": "number"},
"GPT_Summary": {"type": "string"}
},
"required": [
"Rank",
"URL",
"Probability",
"GPT_Summary"
]
}
}
}
)
return json.loads(response.choices[0].message.content)["results"]

We use the response_format option to ensure the output is returned as a structured JSON object, making parsing and downstream processing much easier. Be aware that the OpenAI API has usage quotas and token limits that depend on your subscription tier and the model in use. Occasionally, the API may experience outages or degraded performance — always check OpenAI’s status page or the official documentation if things go wrong. For this project, we use the gpt-4o model, which offers faster responses and lower cost per token. Keep in mind that token limits and output behavior can vary across different models and versions.

Step 4: Full Pipeline Example

Now let’s combine everything into one pipeline and try it out:

query = "best waterproof trail running shoes"
raw = fetch_google_results(query, serpapi_key)
df = clean_results(raw)
response = query_openai(df, keyword="trail running shoes")
print(pd.DataFrame(response).head())

## Rank ... GPT_Summary
## 0 1 ... Guide on top waterproof trail running shoes.
## 1 2 ... Listing of waterproof trail running shoes at REI.
## 2 3 ... Reddit discussion on waterproof trail running ...
## 3 4 ... Review of best waterproof trail running shoes ...
## 4 5 ... Tested reviews of top trail running shoes for ...
##
## [5 rows x 4 columns]

if response:
df_scores = pd.DataFrame(response)
final_df = pd.merge(df, df_scores, on=["Rank", "URL"],how="inner")
print(final_df[["Rank", "Title", "Probability","GPT_Summary"]].head())
else:
print("GPT analysis failed.")

## Rank ... GPT_Summary
## 0 1 ... Guide on top waterproof trail running shoes.
## 1 2 ... Listing of waterproof trail running shoes at REI.
## 2 3 ... Reddit discussion on waterproof trail running ...
## 3 4 ... Review of best waterproof trail running shoes ...
## 4 5 ... Tested reviews of top trail running shoes for ...
##
## [5 rows x 4 columns]

Im only going to use R for the final table, but you can use any other language to do it.

knitr::kable(
py$final_df[, c("Rank", "URL","Title", "Probability", "GPT_Summary")],
format = "html",
caption = "GPT Analysis of Search Results"
)

Table 1: GPT Analysis of Search Results
Rank	URL	Title	Probability	GPT_Summary
1	https://runrepeat.com/guides/best-trail-waterproof-running-shoes	7 Best Waterproof Trail Running Shoes in 2025	0.95	Guide on top waterproof trail running shoes.
2	https://www.rei.com/c/trail-running-shoes/f/f-waterproof	Waterproof Trail-Running Shoes \| REI Co-op	0.90	Listing of waterproof trail running shoes at REI.
3	https://www.reddit.com/r/trailrunning/comments/1adsza0/need_some_recommendations_for_waterproof_running/	Need some recommendations for waterproof running ...	0.75	Reddit discussion on waterproof trail running shoe recommendations.
4	https://www.livefortheoutdoors.com/trail-running/shoes/best-waterproof-trail-running-shoes/	Best waterproof trail running shoes for 2025 \| Tested and ...	0.95	Review of best waterproof trail running shoes for 2025.
5	https://www.outdoorgearlab.com/topics/shoes-and-boots/best-trail-running-shoes-womens	The 9 Best Trail Running Shoes for Women \| Tested	0.85	Tested reviews of top trail running shoes for women.
6	https://www.brooksrunning.com/en_us/mens/shoes/trail-shoes/?srsltid=AfmBOoow-beP7iwwQgxtcrDarJJkMzB4EGLSRJQnxacCDaeJa8C2rU-a	Men's Trail Running Shoes	0.80	Brooks' men's trail shoes with GORE-TEX options.
7	https://www.rei.com/learn/expert-advice/should-you-get-waterproof-trail-running-shoes.html	Should I Get Waterproof Trail-Running Shoes	0.90	Guide on deciding to purchase waterproof trail running shoes.
8	https://www.trailandkale.com/roundups/best-trail-running-shoes-buyers-guide/	The Best Trail Running Shoes In 2025 [Tried And Tested]	0.95	Comprehensive guide to best trail running shoes.
9	https://singletrackworld.com/forum/off-topic/waterproof-trail-running-shoes-yay-or-nay/	Waterproof trail running shoes – Yay or nay? – Chat Forum	0.70	User opinions on waterproof trail shoes in forum.
10	https://lifehacker.com/health/waterproof-running-shoes-are-game-changer	I'm a Runner, and These Water-Resistant Running Shoes ...	0.80	Review of water-resistant running shoes for runners.

Conclusion

In this post, we built a modular and extensible web search analyzer using Python, SerpAPI, and OpenAI’s GPT. We demonstrated how to automatically retrieve search results, clean and structure the data, and use a large language model to generate useful summaries and rank the most relevant links based on probabilistic reasoning.

This framework is highly adaptable:

You can modify the SerpAPI parameters to target different locations, languages, or result depths.
You can enhance the GPT prompt to fit specific domains — whether it’s shopping, education, or scientific research.
The structured outputs can be stored in CSVs, Excel files, or databases for further visualization or batch processing.
You can even integrate this pipeline into larger workflows, such as product monitoring, competitive analysis, or content recommendation engines.

By abstracting complexity through APIs and combining them with thoughtful logic in Python, we enable powerful automation with very little code. This project is a good example of how language models and APIs can work together to extract meaning and prioritize knowledge at scale — without writing a single web scraper.

Happy exploring!

An Introduction to Regression Discontinuity Design (RDD)

Thu, 06 Feb 2025 00:00:00 +0000

Introduction

Have you ever wondered how a government might evaluate the impact of a new scholarship program that only kicks in when a student’s GPA exceeds a specific threshold? Or how a policy change affecting tax rates for households above a certain income level actually influences economic behavior? These are just a couple of scenarios where Regression Discontinuity Design (RDD) comes into play—a powerful tool for causal inference when treatment assignment hinges on a cutoff.

In this blog post, we’ll explore RDD by building a local linear regression from the ground up—no magic packages or black-box functions here! Instead, we’ll simulate data that features a clear jump at a predetermined cutoff and then walk through the entire estimation process step by step. You’ll see exactly how the treatment effect is estimated right at the threshold and how the confidence interval is constructed to quantify uncertainty.

By starting from zero, we aim to demystify the inner workings of RDD, ensuring that you not only learn how to implement it but also understand the underlying logic behind the method. Whether you’re evaluating educational policies, economic interventions, or any scenario with a clear threshold effect, this hands-on guide is your first step toward mastering RDD.

Regression Discontinuity Design (RDD) is a quasi-experimental method used to estimate causal effects when treatment assignment is determined by a cutoff in a continuous variable (the running variable). The key idea is that units just below and just above the cutoff are very similar, so a jump in the outcome at the cutoff can be interpreted as the causal effect of the treatment.

Mathematically, let $$ m(x) = E[Y \mid X=x] $$

be the conditional expectation of the outcome Y given the running variable X. Suppose that there is a discontinuity at x = c so that: $$ m(x) = \begin{cases} m_0(x), & x < c, \newline m_1(x), & x \ge c. \end{cases} $$ The treatment effect at the cutoff is then estimated as: $$ \hat{\tau} = m_1(c) - m_0(c). $$

The variance of the estimated effect is: $$ \operatorname{Var}(\hat{\tau}) = \operatorname{Var}(m_1(c)) + \operatorname{Var}(m_0(c)) $$

and the 95% confidence interval (CI) is given by: $$ \hat{\tau} \pm z_{0.025} \sqrt{\operatorname{Var}(m_1(c)) + \operatorname{Var}(m_0(c))} $$ with $ z_{0.025} \approx 1.96$.

In the following sections, we illustrate these concepts with a simulated dataset.

1. Simulated Data

We simulate a dataset with 500 observations. The running variable x is drawn uniformly from [-10,10]. The outcome y is generated with a discontinuity (a jump) at x=0. In our simulation, for x < 0 we set: $$ y = 5 + 0.3x + \varepsilon, $$ and for $x \ge 0$: $$ y = 7 + 0.3x + \varepsilon, $$ with$\varepsilon \sim N(0,1)$.

import numpy as np
import pandas as pd
# Set the seed for reproducibility
np.random.seed(42)
# Simulate 500 observations for the running variable x
N = 500
x = np.random.uniform(-10, 10, size=N)
# Generate noise
error = np.random.normal(0, 1, size=N)
# Generate outcome y with a jump at x = 0 (jump of 2 units)
y = np.where(x < 0, 5 + 0.3 * x, 7 + 0.3 * x) + error
# Create a DataFrame
data = pd.DataFrame({'x': x, 'y': y})
print(data.head())

## x y
## 0 -2.509198 4.588997
## 1 9.014286 11.580457
## 2 4.639879 9.342387
## 3 1.973170 7.015047
## 4 -6.879627 2.037697

In this simulated example, the cutoff is at $x=0$.

2. Epanechnikov Kernel

The Epanechnikov kernel is used to assign weights to observations based on their distance from the cutoff. Observations closer to the cutoff are given higher weights. The kernel is defined as: $$ K(u) = \begin{cases} 0.75 (1 - u^2) & \text{if } |u| \le 1, \newline 0 & \text{otherwise.} \end{cases} $$ This topic is explained more deeply on my other blog post about Kernel Density Estimation.

def epanechnikov_kernel(u):
"""
Epanechnikov kernel function.
Returns:
K(u) = 0.75 * (1 - u^2) for |u| <= 1, and 0 otherwise.
"""
if abs(u) <= 1:
return 0.75 * (1 - u**2)
else:
return 0

3. Local Linear Regression

To estimate the outcome function on either side of the cutoff, we perform a weighted linear regression. We center the running variable at the cutoff $c=0$ so that: $$ x_{\text{centered}} = x - c = x. $$ The local linear regression model is: $$ y = \beta_0 + \beta_1 x_{\text{centered}} + \varepsilon. $$ At $x_{\text{centered}} = 0$, the predicted outcome is $\beta_0$, which estimates m(c).

def local_linear_regression(df, y_variable, x_variable, cutoff, bandwidth, kernel_function):
"""
Performs local linear regression on observations below and above the cutoff.
The running variable is centered at the cutoff:
$$ x_{centered} = x - c $$
so that the intercept represents the estimated outcome at the cutoff, i.e., $$ m(c) = \\beta_0 $$.
"""
df_local = df.copy()
# Center the running variable at the cutoff (c = 0)
df_local['x_centered'] = df_local[x_variable] - cutoff
df_local['constant'] = 1
# Compute kernel weights based on the scaled distance from the cutoff
df_local['weights'] = df_local['x_centered'].apply(lambda x: kernel_function(x / bandwidth))
# Split data into observations below and above the cutoff
df_below = df_local[df_local['x_centered'] <= 0]
df_above = df_local[df_local['x_centered'] > 0]
def _weighted_regression(sub_df, y_var):
if sub_df.empty:
return np.array([np.nan, np.nan]), np.full((2, 2), np.nan)
X = sub_df[['constant', 'x_centered']].values
Y = sub_df[y_var].values
W = np.diag(sub_df['weights'].values)
try:
XTWX_inv = np.linalg.inv(X.T @ W @ X)
theta = XTWX_inv @ (X.T @ W @ Y)
residuals = Y - X @ theta
variance_matrix = XTWX_inv @ (X.T @ W @ np.diag(residuals**2) @ W @ X) @ XTWX_inv
except np.linalg.LinAlgError:
return np.array([np.nan, np.nan]), np.full((2, 2), np.nan)
return theta, variance_matrix
theta_below, var_below = _weighted_regression(df_below, y_variable)
theta_above, var_above = _weighted_regression(df_above, y_variable)
return theta_below, var_below, theta_above, var_above

RDD Estimator

The RDD estimator computes the treatment effect as the difference in the predicted outcomes at the cutoff: $$ \hat{\tau} = \hat{m}_1(0) - \hat{m}_0(0), $$ where $\hat{m}_1(0)$ and $\hat{m}_0(0)$ are the intercepts from the regressions above and below the cutoff, respectively. The variance of the treatment effect is: $$ \operatorname{Var}(\hat{\tau}) = \operatorname{Var}(\hat{m}_1(0)) + \operatorname{Var}(\hat{m}_0(0)), $$ and its standard error is the square root of this sum.

def rdd_estimator(df, y_variable, x_variable, cutoff, bandwidth, kernel_function):
"""
Estimates the RDD treatment effect as the difference in the predicted outcome at the cutoff.
Since the running variable is centered at the cutoff,
the predicted outcome is simply the intercept:
$$ \\hat{m}(c) = \\beta_0. $$
Thus, the treatment effect is:
$$ \\hat{\\tau} = \\hat{m}_1(c) - \\hat{m}_0(c). $$
"""
theta_below, var_below, theta_above, var_above = local_linear_regression(
df, y_variable, x_variable, cutoff, bandwidth, kernel_function)
if np.any(np.isnan(theta_below)) or np.any(np.isnan(theta_above)):
return np.nan, np.nan
# Predicted outcome at the cutoff is the intercept (when x_centered = 0)
pred_below = theta_below[0]
pred_above = theta_above[0]
treatment_effect = pred_above - pred_below
# Variance of the intercept estimates (using the selector vector e0 = [1, 0])
e0 = np.array([1.0, 0.0])
se = np.sqrt(e0 @ (var_below + var_above) @ e0)
return treatment_effect, se

Running the Analysis and Constructing the Confidence Interval

We now run the RDD estimator using a chosen bandwidth and compute the 95% confidence interval. Recall, the 95% CI is given by: $$ CI_\text{95%} = \hat{\tau} \pm 1.96 \times SE(\hat{\tau}). $$

from scipy.stats import norm
# Choose a bandwidth (experiment with this value)
bandwidth = 2.0
# Run the RDD estimator on the simulated data. Here, the cutoff is 0.
treatment_effect, se = rdd_estimator(data, 'y', 'x', cutoff=0, bandwidth=bandwidth, kernel_function=epanechnikov_kernel)
if not np.isnan(treatment_effect):
print("RDD Estimate (ATE): {:.4f}".format(treatment_effect))
print("Standard Error: {:.4f}".format(se))
# Compute the 95% confidence interval.
alpha = 0.05
z = norm.ppf(1 - alpha/2) # Approximately 1.96 for a 95% CI.
ci_lower = treatment_effect - z * se
ci_upper = treatment_effect + z * se
print("95% Confidence Interval: [{:.4f}, {:.4f}]".format(ci_lower, ci_upper))
else:
print("RDD estimation failed. Please check your data or bandwidth.")

## RDD Estimate (ATE): 1.9539
## Standard Error: 0.3026
## 95% Confidence Interval: [1.3607, 2.5470]

6. Graphical Representation: Plotting the Local Linear Fits and Confidence Intervals

Below we plot the estimated local linear regressions on either side of the cutoff, along with their pointwise 95% confidence intervals.

import matplotlib.pyplot as plt
def plot_rdd_fit(df, y_variable, x_variable, cutoff, bandwidth, kernel_function, n_grid=100, alpha=0.05):
"""
Plots the local linear regression fits with pointwise 95% confidence intervals.
The regression is estimated separately on each side of the cutoff.
"""
# Get the estimates from local linear regressions.
theta_below, var_below, theta_above, var_above = local_linear_regression(df, y_variable, x_variable, cutoff, bandwidth, kernel_function)
if np.any(np.isnan(theta_below)) or np.any(np.isnan(theta_above)):
print("Local-linear regression failed on one side. No plot produced.")
return
# Extract coefficients for each side.
beta0_below, beta1_below = theta_below
beta0_above, beta1_above = theta_above
# Create grids for values below and above the cutoff.
x_min, x_max = df[x_variable].min(), df[x_variable].max()
grid_below = np.linspace(x_min, cutoff, n_grid)
grid_above = np.linspace(cutoff, x_max, n_grid)
# Compute the z-value for the (1-alpha) CI.
z_val = norm.ppf(1 - alpha/2)
# Compute predictions and confidence intervals for observations below the cutoff.
mhat_below = []
ci_lower_below = []
ci_upper_below = []
for xg in grid_below:
# Prediction: m0(x) = beta0_below + beta1_below * (x - cutoff)
pred = beta0_below + beta1_below * (xg - cutoff)
vec = np.array([1.0, xg - cutoff])
var_pred = vec @ var_below @ vec
se_pred = np.sqrt(var_pred)
mhat_below.append(pred)
ci_lower_below.append(pred - z_val * se_pred)
ci_upper_below.append(pred + z_val * se_pred)
# Compute predictions and confidence intervals for observations above the cutoff.
mhat_above = []
ci_lower_above = []
ci_upper_above = []
for xg in grid_above:
# Prediction: m1(x) = beta0_above + beta1_above * (x - cutoff)
pred = beta0_above + beta1_above * (xg - cutoff)
vec = np.array([1.0, xg - cutoff])
var_pred = vec @ var_above @ vec
se_pred = np.sqrt(var_pred)
mhat_above.append(pred)
ci_lower_above.append(pred - z_val * se_pred)
ci_upper_above.append(pred + z_val * se_pred)
#Lets create the original lines from the formula
y_below = 5 + 0.3 * grid_below
y_above = 7 + 0.3 * grid_above
# Plot the results.
plt.figure(figsize=(8,6))
plt.plot(grid_below, mhat_below, color='blue', label=r'$\hat{m}_0(x)$')
plt.fill_between(grid_below, ci_lower_below, ci_upper_below, color='blue', alpha=0.3)
plt.plot(grid_above, mhat_above, color='red', label=r'$\hat{m}_1(x)$')
plt.fill_between(grid_above, ci_lower_above, ci_upper_above, color='red', alpha=0.3)
plt.axvline(cutoff, color='black', linestyle='--', label='Cutoff')
plt.xlabel('Running Variable')
plt.ylabel('Outcome')
plt.title('Local Linear Regression Fit with 95% Confidence Intervals')
#Lets also add the original data
plt.scatter(df[x_variable], df[y_variable], color='black', alpha=0.15, label='Data')
#Lets plot the original lines
plt.plot(grid_below, y_below, color='green', label='y = 5 + 0.3x',linestyle='-.')
plt.plot(grid_above, y_above, color='purple', label='y = 7 + 0.3x',linestyle='-.')
plt.legend()
plt.show()
# Call the plotting function.
plot_rdd_fit(data, 'y', 'x', cutoff=0, bandwidth=bandwidth, kernel_function=epanechnikov_kernel)

As we can see, the local linear regression effectively emphasizes the data closest to the cutoff by using weighted observations. While the slope estimates may be more variable—since they rely on data closer to the cutoff, the intercept provides a robust estimate of the treatment effect at the cutoff. This is why the confidence interval is narrow near the cutoff (where the model is most precise) and widens as we move away, reflecting that our inference is primarily focused on the discontinuity point. Overall, the model accurately captures the differences in outcomes between the two sides of the cutoff.

Conclusion

In this post we have:

Simulated data with a discontinuity at $x=0$,
Used an Epanechnikov kernel to weight observations based on their proximity to the cutoff,
Performed local linear regressions on each side of the cutoff to estimate the outcome functions, and
Estimated the treatment effect as the jump at the cutoff, along with its 95% confidence interval.

This example provides a didactic introduction to RDD. By focusing on data near the cutoff, RDD isolates the causal effect of the treatment. Happy coding and exploring RDD!

Analyzing Bitcoin Options with Deribit API and Black-Scholes Model

Tue, 26 Nov 2024 00:00:00 +0000

Introduction

Bitcoin options have become an essential tool in modern financial markets, offering traders and researchers opportunities to manage risk, speculate on price movements, and explore market dynamics. However, understanding and leveraging these instruments require a combination of financial knowledge and technical expertise.

This blog aims to provide a step-by-step guide to analyzing Bitcoin options by integrating key financial concepts with practical coding implementations. Specifically, we will:

Download and Process Historical Data:
- Learn how to use the Deribit API to retrieve detailed Bitcoin options data, overcoming common challenges in data collection through Python automation.blog
Understand the Black-Scholes Model:
- Explore the mathematical foundations of options pricing, focusing on how variables like volatility and time-to-maturity influence option values.
Apply Machine Learning for Volatility Forecasting:
- Use a Random Forest model to predict implied volatility, showcasing a modern, data-driven approach to financial analysis.
Implement a Delta Strategy:
- Discover how delta, a key measure in options trading, can be used to adjust portfolios dynamically and manage risk effectively.

We will be showing:

API expertise: Using Python to interact with APIs and process large datasets effectively.
Finance Foundations: Applying theoretical models like Black-Scholes in real-world scenarios.
Machine Learning: Leveraging data-driven approaches to forecast critical metrics in financial markets.
Strategic Thinking: Translating insights from data into actionable strategies, like delta crypto finance strategies

Part 1: Downloading Historical Bitcoin Options Data

We will use Deribit’s historical API to access Bitcoin options data. This API allows fetching complete historical data for options trades.

The standard Deribit API imposes a limitation on accessing historical data, allowing users to retrieve trades only up to two weeks in the past. To overcome this, Deribit provides a separate historical API endpoint, enabling access to older trade data with detailed records of options activity.

The steps involve setting up a loop to retrieve data in manageable time intervals and processing the results for further analysis.

Code Implementation:

import requests
import pandas as pd
from datetime import datetime, timedelta
import time
# Define the historical trades URL
historical_trades_url = "https://history.deribit.com/api/v2/public/get_last_trades_by_currency_and_time"
# Function to get historical trades
def get_historical_trades(currency, kind, start_time, end_time, count=1000):
# Convert timestamps to milliseconds since epoch
end_timestamp = int(time.mktime(end_time.timetuple()) * 1000)
start_timestamp = int(time.mktime(start_time.timetuple()) * 1000)
# Define the payload for fetching historical trades
trades_payload = {
'currency': currency,
'start_timestamp': start_timestamp,
'end_timestamp': end_timestamp,
'kind': kind,
'count': count,
'include_old': True
}
# Send the request to get historical trades
trades_response = requests.get(historical_trades_url, params=trades_payload)
# Check if the response is successful
if trades_response.status_code == 200:
trades_data = trades_response.json()
#print("Historical Trades Response:", trades_data)
# Process and print trades data
if 'result' in trades_data:
print("Successfully fetched trades data")
else:
print("Failed to fetch trades data")
else:
print(f"Failed to fetch data: {trades_response.status_code}")
return trades_data
# Function to process and filter trades data
def get_trades(trades_data):
trades = trades_data['result']['trades']
df = pd.json_normalize(trades)
if 'timestamp' in df.columns:
df['timestamp'] = pd.to_datetime(df['timestamp'], unit='ms')
# Filter for call options
calls = df[df['instrument_name'].str.endswith('C')].copy()
# Extract maturity and strike prices
calls.loc[:, 'maturity'] = calls['instrument_name'].apply(lambda x: x.split('-')[1])
calls.loc[:, 'strike'] = calls['instrument_name'].apply(lambda x: x.split('-')[2].replace('C', '').replace('P', ''))
calls['strike'] = calls['strike'].astype(float)
calls['maturity'] = pd.to_datetime(calls['maturity'], format='%d%b%y')
else:
# print("No timestamp column for date", start_time)
calls = pd.DataFrame()
return calls
# Loop to fetch data for multiple time intervals
historical_trades = []
start_time = datetime(2024, 11, 7) # Define starting time
for i in range(12): # Loop to get data for 24 hours
end_time = start_time + timedelta(hours=1)
trades_response = get_historical_trades('BTC', 'option', start_time, end_time)
historical_trades.append(get_trades(trades_response))
start_time = end_time
time.sleep(3) # Wait 3 seconds between API calls

## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data
## Successfully fetched trades data


historical_trades = pd.concat(historical_trades, ignore_index=True)
print(historical_trades.head())

## trade_seq trade_id ... combo_trade_id liquidation
## 0 679 324402707 ... NaN NaN
## 1 43 324402703 ... NaN NaN
## 2 30 324402702 ... NaN NaN
## 3 678 324402688 ... NaN NaN
## 4 4318 324402676 ... NaN NaN
##
## [5 rows x 18 columns]

Here’s a detailed explanation of what happens step by step:

Connecting to Deribit’s Historical API

The script interacts with the Deribit historical API endpoint, which allows fetching trade data for Bitcoin options. To use this API effectively, the script sends a structured request containing details such as the time range for the data, the type of instrument (e.g., options), and the currency (Bitcoin).

Data Retrieval in Time Intervals

Since the historical database can be extensive, the script divides the data request into manageable chunks. Specifically, it retrieves data in 1-hour intervals, starting from a defined start_time. This approach ensures that the API isn’t overwhelmed by large requests and adheres to any rate limits imposed by the server. After each request, the script waits for a few seconds before proceeding to the next, ensuring smooth communication with the API.

Filtering and Processing the Data

Once the trade data is retrieved, it is processed to extract relevant information. The script focuses on “call options,” a specific type of option contract, by filtering instruments whose names end with “C.” It further extracts important details such as:

Strike Price: The agreed price at which the option can be exercised.
Maturity Date: The expiration date of the option.

Automating the Workflow

To handle larger datasets, the script automates the entire workflow using a loop. For each iteration, it fetches data for a 1-hour window, processes the retrieved trades, and appends the results to a list. This list represents a cumulative collection of all the data retrieved across the defined time period (e.g., 7 hours in this example).

In summary:

Historical API: The URL https://history.deribit.com/api/v2/public/get_last_trades_by_currency_and_time allows accessing historical options data. Unlike the WebSocket API, this endpoint provides complete historical data.
Data Retrieval Loop: To avoid overwhelming the server, we query data in 1-hour intervals with a delay between requests.
Filtering and Processing:
- Filters trades to extract only call options (instrument_name ends with C).
- Extracts additional details like strike prices and maturity dates for analysis.

This code ensures efficient fetching of historical Bitcoin options data while adhering to API rate limits and provides well-structured data ready for modeling and analysis.

Part 2: Understanding the Black-Scholes Model

The Black-Scholes model is a mathematical formula for pricing European call and put options. It assumes constant volatility and no dividends.

Essentially, it determines the fair value of a call or put option based on several key factors: the current price of the underlying asset ($S_0$), the option’s strike price ($K$), the time left until expiration ($T$), the risk-free interest rate ($r$), and the volatility of the asset’s price ($ \sigma $). These variables are combined into a formula that accounts for the probabilistic nature of market movements and the time value of money.

At its core, the model calculates the likelihood of an option finishing “in the money” (profitable) by using two probability factors, $N(d_1)$ and $N(d_2)$, derived from the normal distribution. These terms represent the probabilities of different market scenarios. For a call option, the formula calculates the expected payoff by weighing the probability of the stock exceeding the strike price at expiration against the present value of the strike price itself. This balance ensures that the price reflects both the potential reward and the risks involved.

One of the key strengths of the Black-Scholes model is its simplicity and assumptions, such as constant volatility and no dividends. While these assumptions don’t always perfectly match reality, the model is an essential starting point for understanding option pricing and has paved the way for more advanced variations used in modern finance.

Formula:

The price of a call option is given by: $$ C = S_0 \cdot N(d_1) - K \cdot e^{-rT} \cdot N(d_2) $$ Where:

$ S_0 $: Current stock price
$ K $: Strike price
$ T $: Time to expiration
$ r $: Risk-free interest rate
$ N(d) $: Cumulative standard normal distribution
$ d_1 $: $ \frac{\ln(S_0 / K) + (r + \sigma^2 / 2)T}{\sigma \sqrt{T}} $
$ d_2 $: $d_1 - \sigma \sqrt{T}$

Code:

from scipy.stats import norm
import numpy as np
# Black-Scholes Formula
def black_scholes(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
elif option_type == "put":
return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
# Example usage
S = 20000 # Current BTC price
K = 21000 # Strike price
T = 30 / 365 # Time to expiration in years
r = 0.01 # Risk-free rate
sigma = 0.5 # Volatility
call_price = black_scholes(S, K, T, r, sigma)
print(f"Call Option Price: {call_price}, for Strike Price: {K} and Time to Expiration: {T} years")

## Call Option Price: 744.2789814617408, for Strike Price: 21000 and Time to Expiration: 0.0821917808219178 years

Implied Volatility

Implied volatility (IV) is a measure of the market’s expectation of the future volatility of the underlying asset’s price. Unlike historical volatility, which is based on past price movements, IV is derived from the market price of options using models like Black-Scholes. It reflects the degree of expected price fluctuations of the underlying asset over the life of the option.

In practice, implied volatility is an essential tool for traders because:

Market Sentiment: It captures the market’s expectations and risk perceptions.
Option Pricing: Higher IV leads to higher option premiums since the probability of significant price movement increases.
Volatility Smile: In real markets, IV often varies with the option’s strike price, forming a “smile” shape when plotted, contrary to Black-Scholes’ assumption of constant IV.

A “volatility smile” occurs because deep in-the-money and out-of-the-money options often have higher implied volatilities than at-the-money options. This effect reflects market imperfections, demand for certain options, and potential risk.

Black-Scholes Model and Implied Volatility

The price of a call option under the Black-Scholes formula is given by:

$$ C = S_0 \cdot N(d_1) - K \cdot e^{-rT} \cdot N(d_2) $$

Here, $ \sigma $ (volatility) is typically an input into the model to calculate the option’s price. However, in the case of implied volatility, we reverse the process: we use the observed market price of the option to solve for $ \sigma $. This involves numerical optimization, as there is no closed-form solution for $\sigma $.

Limitations of Black-Scholes Implied Volatility

Constant Volatility Assumption:
- The Black-Scholes model assumes that volatility remains constant over the life of the option. However, in reality, volatility is dynamic and varies depending on factors like market conditions, economic data releases, and investor sentiment.
Lognormal Price Distribution:
- The model assumes that asset prices follow a lognormal distribution, meaning they have a symmetric behavior. However, real-world markets often exhibit skewness (asymmetric returns) and kurtosis (fat tails).
Short-Term Maturity Issues:
- For options close to expiration, small changes in price or volatility can lead to large variations in implied volatility. This amplifies numerical optimization errors.
Deep ITM or OTM Options:
- Deep in-the-money (ITM) or out-of-the-money (OTM) options often result in extreme or nonsensical implied volatilities because their market prices are less sensitive to volatility.

Why Does the Market Show a Volatility Smile?

Market Skewness and Risk Aversion:
- Investors often seek protection against extreme moves in the market, leading to higher demand (and hence, higher implied volatility) for deep OTM options, especially put options. This demand reflects their willingness to pay more for tail risk hedging.
Supply-Demand Dynamics:
- Market makers adjust the prices of options based on order flow. If there is higher demand for certain strikes, it increases their implied volatility.
Real-World Volatility is Stochastic:
- Unlike the constant volatility assumed by Black-Scholes, real-world volatility is stochastic (i.e., it changes over time) and often exhibits clustering. This adds curvature to the implied volatility plot.
Market Sentiment and Imbalances:
- In bearish markets, implied volatility tends to increase, particularly for OTM puts, as investors rush to buy downside protection. Conversely, in bullish markets, OTM call options might show elevated implied volatilities due to speculative demand.

The smile in the implied volatility plot reflects the limitations of the Black-Scholes model and the realities of market behavior. While Black-Scholes assumes a flat volatility curve, the observed smile arises from demand for risk hedging, market maker adjustments, and the dynamic nature of volatility. Understanding these nuances allows traders and quants to better interpret market signals and deploy strategies like delta hedging, which benefit from accurate volatility forecasting.

Below is the Python code to calculate and plot the implied volatility smile using the downloaded options data.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
df = historical_trades
df['days_to_maturity'] = (df['maturity'] - df['timestamp']).dt.days
# Step 1: Clean and filter the data
# Remove rows with invalid or unrealistic IV values
filtered_calls = df[(df['iv'] > 0) & (df['iv'] < 500)] # Filter IV column for valid values
# Ensure days to maturity is positive
filtered_calls = filtered_calls[filtered_calls['days_to_maturity'] > 0]
# Step 2: Calculate moneyness (index_price / strike)
filtered_calls['moneyness'] = filtered_calls['index_price'] / filtered_calls['strike']
# Step 3: Filter by similar days to maturity (e.g., 7 to 30 days)
similar_maturity_calls = filtered_calls[
(filtered_calls['days_to_maturity'] >= 7) &
(filtered_calls['days_to_maturity'] <= 30)
]
# Step 4: Plot the implied volatility smile using moneyness
plt.figure(figsize=(10, 6))
plt.scatter(similar_maturity_calls['moneyness'], similar_maturity_calls['iv'], alpha=0.7, label='Implied Volatility')
plt.title("Implied Volatility Smile (7-30 Days to Maturity, Moneyness-Based)")
plt.xlabel("Moneyness (Index Price / Strike Price)")
plt.ylabel("Implied Volatility")
plt.axvline(1, color='red', linestyle='--', label='At-the-Money (ATM)')
plt.grid(True)
plt.legend()
plt.show()

The implied volatility provided by Deribit is derived directly from live market data and reflects the collective sentiment, supply-demand dynamics, and real-time risk perceptions of market participants, without the fridge assumptions of the Black Scholes Formula. This makes it a highly reliable and actionable measure, as it inherently accounts for factors that models like Black-Scholes cannot capture, such as market skewness, kurtosis, and stochastic volatility.

By using Deribit’s IV, we can focus on applying this market-driven insight to strategies like delta hedging, rather than spending time troubleshooting computational issues or compensating for the limitations of traditional models.

Model Deribit Implied Volatility

Lets try to create a model that would forecast the implied volatility of Deribit, using as inputs the characteristics of the option. This model should be non parametric as we dont know what specific model what used by Deribit.

To use the Random Forest model for predicting implied volatility (IV), replacing the analytical Black-Scholes-based calculation, we follow these steps:

Filter and Clean the Data:
- We remove outlier trades that might distort IV (e.g., combo trades, block trades, liquidation events).
- Clean and engineer features such as moneyness, days to maturity, historical volatility, and moving averages.
- Ensure no missing or invalid values remain in the features or target (iv).
Train a Random Forest Model:
- We use relevant features that are correlated with IV, such as days_to_maturity, moneyness, and volatility indicators.
- Train the Random Forest model to predict the IV (target = 'iv') using a training set.
Predict IV for the Dataset:
- We use the trained model to predict IV for the entire dataset or a test set, then replace the calculated IVs in the dataset with the predicted values for further analysis.
Visualize the Volatility Smile:
- We plot the predicted IVs against the strike prices to observe the smile pattern.

Step-by-Step Code

1. Data Preparation

Start by ensuring your data is clean and ready for modeling.

# Filter trades with valid IV and prices
df = historical_trades.copy()
df['maturity'] = pd.to_datetime(df['maturity'], errors='coerce')
df['timestamp'] = pd.to_datetime(df['timestamp'], errors='coerce')
# Drop rows with invalid dates
df = df.dropna(subset=['maturity', 'timestamp'])
# Calculate days to maturity
df['days_to_maturity'] = (df['maturity'] - df['timestamp']).dt.days
df['time_to_maturity'] = df['maturity'] - df['timestamp']
# Filter out invalid IV values and zero index prices
df = df[(df['iv'] > 0) & (df['iv'] < 500)]
df = df[df['index_price'] > 0]
# Create the 'moneyness' column
df['moneyness'] = df['index_price'] / df['strike']
# Extract the hour of the trade
df['hour'] = df['timestamp'].dt.hour
# Ensure non-negative days to maturity
df['days_to_maturity'] = df['days_to_maturity'].apply(lambda x: max(0, x))
# Sort by timestamp and set as index for rolling calculations
df = df.sort_values('timestamp')
df = df.set_index('timestamp')
# Calculate historical volatility and moving averages
df['historical_volatility'] = df['index_price'].rolling('1h').std() / df['index_price'].rolling(window='1h').mean()
df['moving_average'] = df['index_price'] / df['index_price'].rolling('1h').mean()
df['historical_volatility_7'] = df['index_price'].rolling('7h').std() / df['index_price'].rolling(window='7h').mean()
df['moving_average_7'] = df['index_price'] / df['index_price'].rolling('7h').mean()
# Calculate 3-hour rolling metrics
df['historical_volatility_3'] = df['index_price'].rolling(window='3h').std() / df['index_price'].rolling('3h').mean()
df['moving_average_3'] = df['index_price'] / df['index_price'].rolling('3h').mean()
# Reset the index after rolling calculations
df = df.reset_index()
# Drop rows with NaN values after rolling
df = df.dropna(subset=['historical_volatility', 'moving_average', 'historical_volatility_3', 'moving_average_3'])
# Engineer additional features: minutes to maturity and time to maturity in years
df['minutes_to_maturity'] = df['time_to_maturity'].dt.total_seconds() / 60
df = df[df['minutes_to_maturity'] > 0] # Remove rows with non-positive maturity
df['time_to_maturity_year'] = df['minutes_to_maturity'] / (365 * 24 * 60)
# Verify final DataFrame
print(f"Data size after all calculations: {len(df)}")

## Data size after all calculations: 4772

print(df[['moneyness', 'historical_volatility', 'historical_volatility_3', 'moving_average', 'moving_average_3']].head())

## moneyness historical_volatility ... moving_average moving_average_3
## 1 0.934902 0.000415 ... 1.000293 1.000293
## 2 0.997276 0.000353 ... 1.000227 1.000227
## 3 0.959066 0.000345 ... 1.000285 1.000285
## 4 0.959066 0.000325 ... 1.000228 1.000228
## 5 0.831048 0.000292 ... 1.000048 1.000048
##
## [5 rows x 5 columns]

2. Train the Random Forest Model

Select features and train a Random Forest model.

from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
# Select features and target variable
features = df[['days_to_maturity', 'moneyness', 'historical_volatility', 'hour',
'moving_average', 'moving_average_3', 'moving_average_7',
'historical_volatility_3', 'historical_volatility_7']]
target = df['iv']
# Split data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(features, target, test_size=0.2, random_state=42)
# Train the Random Forest model
model = RandomForestRegressor(random_state=42)
model.fit(X_train, y_train)

RandomForestRegressor(random_state=42)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# Predict and evaluate
y_pred = model.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
print(f'Mean Squared Error: {mse}')

## Mean Squared Error: 11.621701853968593

3. Predict IV for the Full Dataset

Use the trained model to predict IV for the entire dataset.

# Predict implied volatility for the entire dataset
df['predicted_iv'] = model.predict(features)

4. Visualize the Implied Volatility Smile

Create a plot of predicted IVs against strike prices.

import matplotlib.pyplot as plt
# Filter dataset for options with 7 to 30 days to maturity
filtered_df = df[(df['days_to_maturity'] >= 7) & (df['days_to_maturity'] <= 30)].dropna(subset=['iv', 'moneyness', 'historical_volatility'])
# Plot IV against moneyness
plt.figure(figsize=(12, 6))
plt.scatter(filtered_df['moneyness'], filtered_df['iv'], color='blue', alpha=0.6, label='Real IV (Deribit)', s=20)
plt.scatter(filtered_df['moneyness'], filtered_df['predicted_iv'], color='orange', alpha=0.6, label='Calculated IV (Random Forest)', s=20)
# Add plot titles and labels
plt.title("Implied Volatility vs. Moneyness (7-30 Days to Maturity)", fontsize=14)
plt.xlabel("Moneyness (Index Price / Strike Price)", fontsize=12)
plt.ylabel("Implied Volatility", fontsize=12)
plt.grid(True)
plt.legend()
# Show the plot
plt.show()

5. (Optional) Feature Importance

Analyze which features are most important in predicting IV.

import numpy as np
import matplotlib.pyplot as plt
# Get feature importance
importances = model.feature_importances_
feature_names = features.columns
# Sort features by importance
sorted_indices = np.argsort(importances)[::-1] # Sort in descending order
sorted_feature_names = feature_names[sorted_indices]
sorted_importances = importances[sorted_indices]
# Plot feature importance
plt.figure(figsize=(10, 6))
plt.barh(sorted_feature_names, sorted_importances, color='skyblue')
plt.xlabel("Importance")
plt.ylabel("Feature")
plt.title("Feature Importance in Predicting IV")
plt.gca().invert_yaxis() # Invert y-axis to show the highest importance at the top
plt.show()

Key Points

Why Random Forest?: It is a fast, flexible and powerful model that can capture non-linear relationships between features and the target (IV). This is useful since implied volatility may not follow a simple analytical formula.
Feature Selection: Ensure that the selected features are relevant and clean to avoid introducing noise into the model.
Interpretability: Use feature importance to understand which variables are driving the model’s predictions.

This approach replaces the analytical Black-Scholes calculation with a data-driven machine learning model to estimate implied volatility. With a bigger data, the random forest could be better trained and more complex models could be made.

Part 3: Implementing a Delta Strategy

What is Delta in Options Trading?

Delta is a measure of how sensitive an option’s price is to changes in the price of the underlying asset. Specifically, for call options, delta represents the expected change in the option’s price for a \1 change in the underlying asset’s price, assuming all other factors remain constant. Mathematically, delta ($\Delta$) is expressed as:

$$ \Delta = \frac{\partial C}{\partial S} $$

Where:

$C$: Option price
$S$: Price of the underlying asset

For call options, delta values range between 0 and 1. A delta of 0.5, for instance, means the option’s price is expected to change by \0.50 for every \1 change in the asset’s price. In the case of put options, delta values range between -1 and 0.

Delta also provides insight into the likelihood of the option expiring in the money (ITM). For example, a call option with a delta of 0.7 suggests a 70% chance of the option being ITM at expiration. Delta evolves dynamically based on factors like the underlying asset’s price, time to maturity, and volatility.

Understanding delta is crucial for constructing hedging strategies, as it allows traders to adjust their positions to manage risk effectively.

Code:

# Delta calculation
def calculate_delta(S, K, T, r, sigma):
d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
return norm.cdf(d1)
# Example: Hedging with 1 option
delta = calculate_delta(S, K, T, r, sigma)
hedge_position = -delta # Short this fraction of BTC
print(f"Delta Position: {hedge_position} BTC")

## Delta Position: -0.39629035723658007 BTC

Understanding Delta and Using it as a Market Signal

Delta evolves dynamically based on the option’s moneyness, time to maturity, and the volatility of the underlying asset. For example:

At-the-Money (ATM): Delta is approximately 0.5.
Out-of-the-Money (OTM): Delta approaches 0 as the option becomes less sensitive to price changes.
In-the-Money (ITM): Delta approaches 1, reflecting greater sensitivity.

By observing changes in delta, traders can infer shifts in the market sentiment and adjust their positions to maximize profits. For instance, when delta increases significantly, it may signal bullish momentum, while a sharp decrease may indicate bearish trends. Delta-based signals can also help traders decide when to buy or sell Bitcoin to align with market conditions effectively.

Strategy: Using Bitcoin Call Options Delta as a Market Signal

One potential strategy involves monitoring the delta of Bitcoin call options to generate buy/sell signals for Bitcoin itself. Here’s how it works:

Set-Up:
- Use a predictive model (e.g., Random Forest) to estimate implied volatility (IV).
- Calculate delta for a specific call option, ensuring it is near the current Bitcoin price (ATM) with moderate time to maturity (e.g., 7–30 days).
Signal Generation:
- Buy Signal: When delta increases significantly ($ $ approaches 1), it may indicate bullish momentum in Bitcoin or future bullish expectation. Buy/sell Bitcoin to benefit from potential upward price movement.
- Sell Signal: When delta decreases sharply ($ $ approaches 0), it may indicate bearish momentum or a future expected fall in price. Sell/buy Bitcoin to avoid potential losses.
Our portfolio value would look like this:

$$ P_{t+1} = 100 + D_{t}(S_{t+1} - S_t) + \dots + D_1(S_2 - S_1) $$

Where:

Initial Cash: The portfolio starts with \100.
Delta ( $D_t $): Represents the sensitivity of the option’s price to changes in the underlying asset’s price at time t.
Stock Price Change ($ S_{t+1} - S_t $): The difference in the underlying asset’s price between time t and t+1.

Execution:
- Adjust the trading position based on delta movements while considering transaction costs and rebalancing thresholds to minimize unnecessary trades.

At every round t, we adjust the portfolio value by buying/selling $D_t * S_{t}$, which would modify the value of our portfolio at t+1 to $ D_t * (S_{t+1} - S_t)$ . This is the core of our strategy, and its succes would depend of updating the value of $ D_t $ at every round t, which at the same time, would be a function of the forecasted implied volatility, as well as the strike price and time to maturity.

Ensuring that $ D_t $ accurately reflects the option’s sensitivity at each time step is crucial.

While the standard way of thinking about options is using the underlying asset’s price $ S_{t}$, for the execution of this strategy, its not the price of the asset that matters, but the delta of the option.

Considering this, $ S_{t}$ in our strategy could be the initial amount of cash of our portfolio. In the first round, as we buy Bitcoin with the whole amount of cash, the value of our portfolio remains unchanged.

Risk Management:
- Use stop-loss orders to limit losses in case of sharp adverse movements.
- Combine delta signals with other indicators (e.g., volatility or moving averages) to enhance decision-making.

This approach provides a practical and dynamic way to leverage options pricing as a tool for directional trading in Bitcoin markets, making delta a valuable market signal for maximizing profitability.

Workflow for Delta Strategy with Random Forest IV

Delta strategy Logic:
- Adjust the stock position dynamically to match the delta of the option (IV forecasted using the Random Forest model).
- Track portfolio value, cash position, and transaction costs over time.
Simulate Portfolio Performance:
- Simulate the performance of the portfolio as you hedge delta over time.
Plot Results:
- Visualize key metrics like the portfolio value over time or transaction costs.

df.iterrows()

## <generator object DataFrame.iterrows at 0x15ac4fab0>

Code

import pandas as pd
import matplotlib.pyplot as plt
# Define the delta hedging function
def delta_hedging(df, forecast_iv_column, initial_stock_price, r=0.01, trading_fee_rate=0.002):
# Initialize portfolio and stock positions
portfolio_value = 0
cash_position = initial_stock_price # Starting with cash equal to initial stock price
stock_position = 0 # Starting with no stock initially
results = []
for index, row in df.iterrows():
S = row['index_price']
K = row['strike']
T = row['days_to_maturity'] / 365 # Convert days to years
date = row['date']
#print(K)
sigma_forecast = row[forecast_iv_column] / 100 # Convert percentage to decimal
# Calculate forecasted delta
delta_forecast = (calculate_delta(S, K, T, r, sigma_forecast))*(-1) #is negative because we are shorting the call
# Calculate the number of stocks to buy/sell to adjust to the new delta
delta_change = delta_forecast - stock_position
transaction_cost = 0
if delta_change != 0:
# Calculate transaction cost: trading fee + network fee
trading_fee = abs(delta_change) * S * trading_fee_rate
#transaction_cost = trading_fee + network_fee*S
transaction_cost = 0
if abs(delta_change * S) > transaction_cost:
# Update cash position considering transaction cost
cash_position -= transaction_cost #we substract the transaction cost from the cash position
# Adjust stock position to maintain delta neutrality
cash_position += delta_change * S #if delta_change is positive, we are buying stock, if negative we are selling stock.
#If we buy stock, we substract the stock price from the cash position
stock_position = delta_forecast #we set the stock position to the forecasted delta
# Calculate actual option price using Black-Scholes
#option_price_actual = black_scholes_call(S, K, T, r, sigma_actual)
# Update portfolio value
portfolio_value = (stock_position * S*(-1)) + cash_position #portfolio value is the stock position times the stock price plus the cash position
# Store results
results.append({
'timestamp': row['timestamp'], # Ensure timestamp exists
'portfolio_value': portfolio_value,
'index_price': S
})
return pd.DataFrame(results)
# Apply delta hedging strategy to the data
df['date']= df['timestamp'].dt.date
results = delta_hedging(df, 'predicted_iv', initial_stock_price=100)

## <string>:3: RuntimeWarning: divide by zero encountered in scalar divide

# Create a dual-axis plot
fig, ax1 = plt.subplots(figsize=(12, 6))
# Left y-axis: Index Price
ax1.plot(results['timestamp'], results['index_price'], label='Index Price', color='orange')
ax1.set_xlabel("Timestamp", fontsize=12)
ax1.set_ylabel("Index Price", fontsize=12, color='orange')
ax1.tick_params(axis='y', labelcolor='orange')
ax1.grid(True)
# Right y-axis: Portfolio Value
ax2 = ax1.twinx()
ax2.plot(results['timestamp'], results['portfolio_value'], label='Portfolio Value', color='blue')
ax2.set_ylabel("Portfolio Value", fontsize=12, color='blue')
ax2.tick_params(axis='y', labelcolor='blue')
# Title
plt.title("Portfolio Value and Index Price Over Time", fontsize=14)
# Show the plot
fig.tight_layout()
plt.show()

Explanation

Delta Calculation:
- The forecasted delta is calculated using the predicted implied volatility (predicted_iv) from the Random Forest model.
Portfolio Adjustments:
- Adjust stock positions dynamically to match the forecasted delta, considering transaction costs.
Visualization:
- Portfolio Value: Track the total portfolio value over time to assess the effectiveness of the strategy.
- Transaction Costs: Plot cumulative transaction costs to understand the costs
- Delta Changes: Analyze the frequency and magnitude of adjustments required to maintain the delta strategy.

Key Points

Accuracy of IV Forecast: The results will depend heavily on the quality of the forecasted IVs (predicted_iv).
Transaction Costs: Including realistic transaction costs provides a more accurate simulation of the strategy’s performance.
Further Analysis: You can experiment with different time intervals (e.g., rebalancing daily or hourly) to evaluate how often the portfolio should be adjusted.

Delta as a Market Signal and its Impact on the Portfolio

In this approach, the portfolio value is influenced by the deltas of all traded call options, acting as a market signal. Delta behaves like a weighted sum: it is low for out-of-the-money options, high for at-the-money options, and particularly significant for options with high volatility and short days to maturity. Longer-dated options are less impactful due to their lower weight.

Delta dynamically adjusts based on the option’s moneyness and time to maturity. For deep out-of-the-money options, delta approaches 0, limiting losses if the price drops. Conversely, for deep in-the-money options, delta nears 1, amplifying gains if the price rises. Shorter time to maturity further sharpens delta changes, making the portfolio more sensitive to price fluctuations when at-the-money.

While this demonstrates the potential for strategic delta adjustments, the intention of this blog is not to define a trading strategy but to explore delta’s role in portfolio dynamics.

Posible modeling of delta?

It could be possible to model delta, though this would not be followed in this blog, i will write the ideas: Variables that affects $ D_i $:

$$ D_i = \Phi(d_1), \quad \text{where } d_1 = \frac{\ln(S/K) + (r + \sigma^2 / 2)T}{\sigma \sqrt{T}} $$

Where:

Stock Price ($S$): As $S$ increases, $d_1$ increases, and so does $(d_1)$ . A higher $S$ leads to a higher delta for a call option, amplifying the positive impact of $S_{i+1} - S_i$ on the portfolio.
Strike Price ( $K$ ): A lower K (deep in-the-money option) increases $ d_1 $, leading to a higher delta. If the option starts in the money, delta is closer to 1, and the portfolio becomes more sensitive to price changes.
Time to Maturity ( T ): As T decreases, $\sqrt{T}$ decreases, which can increase $d_1 $ (depending on the relationship between $S/K$ and $\sigma $). However, as $T $ , delta converges to 1 if the option is in the money and to 0 if it is out of the money.
Volatility ( $\sigma$ ): Higher volatility increases $ d_1$ when $ S/K > 1 $, leading to a higher delta. In volatile markets, delta can change more significantly, affecting the portfolio positively or negatively depending on the price direction.

Cumulative Sum

$$ \sum_{i=1}^{t} D_i(S_{i+1} - S_i) $$

The overall portfolio value is the sum of all $D_i(S_{i+1} - S_i) $ terms. Each term is influenced by: • Consistency of Positive Movements: If most $ S_{i+1} - S_i > 0$ , and delta $D_i$ is positive, the sum of these terms will increase $P_{t+1}$ . • Magnitude of Delta: The magnitude of $D_i$ amplifies or dampens the effect of $S_{i+1} - S_i$ . Higher delta (closer to 1) leads to greater sensitivity to price changes, increasing the portfolio’s reaction to upward movements.

Considering the implications of normal distribution between the stock price, we could try to use non parametric modeling, like i talked about it this post in the blog.

Conclusion

In this blog, we embarked on a journey through the world of Bitcoin options, using tools and strategies that bridge finance, data science, and machine learning. Let’s summarize what we achieved and the tools that made it possible:

Key Takeaways:

Data Collection with the Deribit API:
- We demonstrated how to use the Deribit API to fetch historical Bitcoin options data efficiently.
- By automating the data retrieval process in Python, we handled the complexities of time-interval-based queries, ensuring well-structured datasets ready for analysis.
Understanding the Black-Scholes Model:
- The blog provided a clear and simplified explanation of the Black-Scholes formula, a cornerstone of options pricing in finance.
- We explored implied volatility, its significance in market sentiment, and limitations of the Black-Scholes model, enhancing the understanding of real-world market dynamics.
Applying Machine Learning for Forecasting:
- By employing a Random Forest model, we showcased a modern approach to predict implied volatility from market data.
- This non-parametric method allowed us to move beyond traditional assumptions and adapt to complex, real-world patterns in options trading.
Delta Strategy for Portfolio Management:
- We implemented a delta strategy, adjusting stock positions dynamically based on the delta of Bitcoin options.
- Through visualization and simulation, we highlighted how delta can act as a market signal and a tool for risk management in financial portfolios.

Electricity Demand in Peru: A Regional Analysis

Mon, 25 Nov 2024 00:00:00 +0000

Motivation

Electricity is a cornerstone of modern life, powering industries, homes, and essential services. Understanding the demand for electricity, particularly at the distribution level, is critical for effective resource planning and optimization. Accurate analysis of electricity demand can help identify usage patterns, forecast future needs, and ensure that energy supply meets demand reliably and sustainably. The code of this project can be seen here.

In this analysis, I focus on the demand of distribution companies (empresas distribuidoras) in Peru, which serve as intermediaries between large power generators and end-users like households and businesses. By studying the demand at the departmental (regional) level, this work provides valuable insights into how energy is consumed across different areas.

The outcomes of such an analysis are highly relevant for:

Policy-making: Supporting government efforts to ensure equitable energy distribution and address regional disparities in electricity access.
Investment decisions: Assisting energy companies in identifying regions that require infrastructure upgrades or additional capacity.
Sustainability: Highlighting opportunities for renewable energy integration and reducing inefficiencies in the power grid.

By combining robust data analysis techniques with domain-specific knowledge, this project demonstrates a practical approach to tackling real-world challenges in the energy sector. This type of analysis is not only adaptable to electricity but can also be extended to other resource management problems, showcasing the versatility of data-driven solutions.

Introduction

Electricity in the country starts being generated at the generators, and then it is transported to both large consumers (such as large mines and factories) and distribution companies (empresas distribuidoras). It is through these distribution companies (empresas distribuidoras) that electricity is “distributed” to regulated customers, such as households and medium or small factories.

The purpose of this analysis is to measure the energy demand of the distribution companies (empresas distribuidoras), organized by department (departamento), in Peru.

The demand of distribution companies (empresas distribuidoras) gives us an indicator of the demand from small and medium-sized businesses, as well as households.

Households consume electricity constantly over time, so variations should be attributed to fluctuations in the demand from small and medium-sized factories/mines/others.

The electrical map (mapa eléctrico) of Peru is available on the COES website.

The distribution of electricity in Peru takes the form of a circulatory system, in which the final branches expand like a tree. This is because the final branches typically use a radial distribution system (sistema de distribución radial).

In the radial system (sistema radial), all the electricity demanded in the area enters through a single point.

For example, in the case of Tumbes, all the lines originate from the Zorritos substation (sub-estación Zorritos) and distribute electricity to the Zarumilla, Tumbes, Zorritos, and Mancora substations (subestaciones Zarumilla, Tumbes, Zorritos y Mancora).

As you might guess, this can lead to the problem of information being recorded more than once. We will explore this issue later.

Source

We will use as data the demand from agents (distributors) published by COES: http://www.coes.org.pe/Portal/DemandaBarras/consulta/index?tipo=2#.

This dataset shows the demand of distribution companies (empresas distribuidoras), measured at each available point (which could be substations, transmission lines, etc.).

The information obtained from COES consists of two parts:

The description of the measurement points, which includes:

i) The 5-digit code of the measurement point.
ii) The distribution company (empresa distribuidora).
iii) The substation/transmission line (subestación/linea de transmisión).
iv) The power or load transported by that point.
v) The geographical location.

Table: Table 1: Header Information (Información de cabecera):

PUNTO.DE.MEDICIÓN	21577	22202	21558	22208
NA	ELECTRO CENTRO	ELECTRO CENTRO	ELECTRO CENTRO	ELECTRO CENTRO
NA	4TP-410	CELDA-22.9KV-AUCAYACU	4TP-202	CELDA - 33KV - EJE MANTARO
NA	22.90	0	69	33
NA	ANDAYCHAGUA ELC	AUCAYACU	AYACUCHO	CAMPO ARMIÑO

The observed data at each point, including the date and time when it was measured. The information is recorded every 30 minutes and is measured in MW. To obtain a unit per hour (MWh), the sum of one hour must be divided by 2.

Table: Table 2: Data measured at each point (Data medida en cada punto):

PUNTO.DE.MEDICIÓN	21577	22202	21558	22208
2020-02-15 00:30:00	0.833	0.96	9.1	0.487
2020-02-15 01:00:00	0.842	0.918	8.61	0.469
2020-02-15 01:30:00	0.771	0.906	8.37	0.463
2020-02-15 02:00:00	0.747	0.89	8.17	0.459
2020-02-15 02:30:00	0.725	0.884	8.15	0.462

Finally, I assign each substation (subestación) to a geographical department (departamento).
This is a personal classification, under which I consider that each distribution company (empresa distribuidora) has an assigned region where it distributes energy. For example, Luz del Sur only distributes energy in Lima.

Other companies (empresas) are assigned to multiple departments (departamentos), so I also use Google Maps and OSINERMING to locate the names of the substations (subestaciones).

Table: Table 3: Organization by Departments (Organización en Departamentos):

PUNTO.DE.MEDICIÓN	DEPARTAMENTO	FECHA.HORA./.SUBESTACIÓN	EMPRESA
21520	ANCASH	PALLASCA	HIDRANDINA
21547	LA LIBERTAD	TAYABAMBA	HIDRANDINA
21644	AREQUIPA	REPARTICIÓN	SEAL
22365	LIMA	SANTA ROSA A.	ENEL DISTRIBUCION PERU
21594	PIURA	CONSTANTE	ELECTRO NOR OESTE

The organization by departments (departamentos) I created can be downloaded here.

One problem with summing this data by department (departamento) is that it can be redundant. As mentioned earlier, some measurement points (puntos de medición) originate from others, meaning the information is recorded two or more times.

A more practical way to observe this issue is by using the single-line diagram (diagrama unifilar) that shows direct connections. This map can also be downloaded from the COES website. Information is only counted more than once in radial systems (sistemas radiales).

On the map, we can see how the demand of the Tumbes, Puerto Pizarro, and Mancora substations (subestaciones) is measured entirely by the line LT-6665A towards the Zorritos - Tumbes 60KV substation (subestación).

We can confirm this in the following table:

## Joining with `by = join_by(PUNTO.DE.MEDICIÓN)`

## Joining with `by = join_by(PUNTO.DE.MEDICIÓN, DEPARTAMENTO,
## `FECHA.HORA./.SUBESTACIÓN`, EMPRESA)`

Table: Table 4: Demand - Tumbes (Demanda - Tumbes)

PUNTO.DE.MEDICIÓN	21588	21585	21586	21589	21590	21614	suma
NA	ELECTRO NOR OESTE	ELECTRO NOR OESTE	ELECTRO NOR OESTE	ELECTRO NOR OESTE	ELECTRO NOR OESTE	ELECTRO NOR OESTE	NA
NA	TR-60/22.9/10kV	TR-60/33/10kV	TR-60/33/10kV	TR-60/22.9/10kV	TR-60/22.9/10kV	LT-6665A	NA
NA	22.90	10	33	10	22.90	60	NA
NA	PUERTO PIZARRO	TUMBES	TUMBES	ZARUMILLA	ZARUMILLA	ZORRITOS - TUMBES	NA
2020-02-15 00:00:00	6.07	9.82	5.62	3.11	1.98	26.6	26.60
2020-02-15 00:30:00	4.86	0.31	5.92	2.97	1.86	15.92	15.92
2020-02-15 01:00:00	4.57	9.1	5.85	2.88	1.77	24.17	24.17
2020-02-15 01:30:00	3.65	9.12	5.83	2.87	1.71	23.18	23.18
2020-02-15 02:00:00	3.68	9.59	5.65	2.57	1.42	22.91	22.91

As seen in both the table and the images above, the measurement point (punto de medición) on the transmission line Zorritos-Tumbes - LT6665A already contains the sum of the points in Tumbes, Puerto Pizarro, and Zarumilla.

Processing

For this presentation, I will analyze two time intervals: 29 days before the start of the coronavirus quarantine (March 16, 2020) and the 29 days following it.

You can download the data I use from the COES website or directly from here.

The information can present a variety of issues:

Missing values (valores omitidos).
Incorrectly entered values (some numbers are multiplied by 10 or more).
Outliers (irrationally different numbers that can persist for an entire day).

Handling missing values

As an example of how to deal with the first type of error (missing values), I will use the data available for Tumbes, according to my categorization, where missing values can be observed.

Visually, it is clear that the series exhibit different seasonality: there are hours during the day when more electricity is consumed and days when less electricity is consumed (such as weekends).

To impute (fill missing values) in the series, we first need to determine the frequency of these seasonal effects and then proceed to impute within each seasonality (for example, this series has a seasonality of 24 hours and another of 7 days). The package ImputeR allows us to do this using the na_seasplit function.

The imputeTS package in R is a powerful tool designed specifically for handling missing values in time series data. One of its key functions, na_seasplit(), is particularly effective for time series with strong seasonal patterns.

In this analysis, I chose na_seasplit() because it imputes missing values by splitting the time series into seasonal components, such as daily or weekly cycles. This approach ensures that the imputed values respect the inherent structure of the data. For example, electricity demand often exhibits predictable patterns based on the time of day or day of the week. By accounting for these patterns, na_seasplit() provides more realistic and accurate imputations compared to generic methods like mean or linear interpolation.

This method is especially useful when:

Missing values occur in a regular seasonal context (e.g., during certain hours or days).
Maintaining the integrity of periodic fluctuations in the data is essential for accurate downstream analysis.

Overall, the use of imputeTS and na_seasplit() demonstrates how specialized tools can enhance the reliability of results, particularly in domains like electricity demand analysis where seasonality plays a crucial role.

The theory behind this process is not the focus of this analysis; however, I recommend the book “Forecasting: Principles and Practice” by Rob Hyndman and his course on Data Camp for further exploration.

As observed in the graph, the missing values (valores omitidos) have been replaced, considering both the daily and weekly seasonality of the series.

A quicker option would be to aggregate the data by day and impute missing values considering only daily seasonality. This will be the method used moving forward, as it helps to avoid errors.

Addressing outliers

To address outliers (valores atípicos), I use the measurements in the Lambayeque department (departamento), according to my classification, as an example. There is an outlier in the measurement point (punto de medición) CAYALTI, which shows an irrational consumption for one day.

Outliers can significantly distort analysis, especially in datasets with numerical values like electricity demand. To avoid this type of error, I use two filters:

i) First, I eliminate any value that exceeds a rational consumption within the 30-minute measurement interval (e.g., 300 MWh).

ii) Second, I eliminate values that exceed 4 times the standard deviation of each series. This approach aims to remove non-recurrent outliers. This method works by identifying data points that deviate excessively from the average pattern of the series, which are likely to be errors or anomalies.

This approach is effective because:

It adapts to the natural variability of each series, ensuring that normal fluctuations are not mistakenly removed.
The 4-standard-deviation threshold strikes a balance between sensitivity and robustness, capturing extreme outliers while retaining most valid data points (on a normal distribution, 99% of the values lie within 3 standart deviations).

By applying this filter, I can clean the data in a systematic way, ensuring that the analysis focuses on meaningful trends without being skewed by extreme and unlikely values.

In the following image, it can be observed how the outliers (valores atípicos) in CAYALTI were removed, and subsequently, the missing value was imputed to fill the series.

However, at the point TUMÁN, a third type of error is still present: for one day, the data was arbitrarily recorded at an incorrect scale (multiplied by 100).

The only way to correct this type of error is to visually inspect the series and manually adjust them.
These types of outliers can affect the interpretation of our data. For example, when analyzing demand by department (departamento), the measurement error in Lambayeque could lead to a misinterpretation.

This code processes electricity consumption data by cleaning, normalizing, and imputing missing values while addressing outliers to prepare it for analysis and visualization. It groups data by date and department, summing electricity demand across all relevant points of measurement for each department. Seasonal imputation (na_seasplit) is applied to ensure smooth trends, and outliers are handled using scaling and thresholding. Finally, it generates a faceted plot showing electricity consumption trends for each department, marking the impact of quarantine (March 15, 2020) and highlighting any anomalies, such as the outlier in Lambayeque. It does this by aggregating all the points of measurement per department.

cont = z %>%
mutate_if(is.character,as.double) %>%
mutate(Fecha = as_date(PUNTO.DE.MEDICIÓN - minutes(30))) %>%
select(-PUNTO.DE.MEDICIÓN) %>%
mutate(Fecha2 = if_else(Fecha > dmy("15-03-2020"),"DC","AC")) %>%
mutate_at(vars(-c(Fecha,Fecha2)),limite) %>%
mutate_at(vars(-c(Fecha,Fecha2)),abs) %>%
group_by(Fecha2) %>%
mutate_at(vars(-c(Fecha,Fecha2)),scale2) %>%
ungroup() %>%
mutate_at(vars(-c(Fecha,Fecha2)),na_ma,maxgap=47)

pordepa= cont %>%
{.[, which(colMeans(!is.na(.)) > 0.5)]} %>%
{.[, which(colMeans(.!= 0,na.rm = TRUE) > 0.5)]} %>%
group_by(Fecha2) %>%
mutate_at(vars(-c(Fecha,Fecha2)), na_seasplit,algorithm="ma") %>%
ungroup() %>%
mutate_at(vars(-c(Fecha,Fecha2)), funs(. / 2000)) %>%
pivot_longer(-c(Fecha,Fecha2),names_to = "PUNTO.DE.MEDICIÓN",values_to = "gwh") %>%
left_join(nombres) %>%
group_by(Fecha,DEPARTAMENTO) %>%
summarise(gwh = sum(gwh)) %>%
ungroup()
factorizar = pordepa %>%
filter(Fecha < as_date("2020-03-15")) %>%
group_by(DEPARTAMENTO) %>%
summarise(gwh = sum(gwh)) %>%
arrange(gwh) %>%
select(DEPARTAMENTO)%>%
tail(7) %>%
{rbind(c("OTROS"),.)}
nombres = mutate(nombres,DEPARTAMENTO = ifelse(nombres$DEPARTAMENTO %in% factorizar[[1]],DEPARTAMENTO,"OTROS"))

cont %>%
{.[, which(colMeans(!is.na(.)) > 0.5)]} %>%
group_by(Fecha2) %>%
ungroup() %>%
mutate_at(vars(-c(Fecha,Fecha2)), funs(. / 2000)) %>%
pivot_longer(-c(Fecha,Fecha2),names_to = "PUNTO.DE.MEDICIÓN",values_to = "gwh") %>%
left_join(nombres) %>%
group_by(Fecha,DEPARTAMENTO) %>%
summarise(gwh = sum(gwh)) %>%
ungroup() %>%
group_by(DEPARTAMENTO) %>%
mutate(gwh = na_seasplit(gwh,find_frequency = TRUE)) %>%
ungroup() %>%
ggplot(aes(x=Fecha, y =gwh)) +
geom_line() +
facet_wrap(.~factor(DEPARTAMENTO,levels = rev(factorizar[[1]])),scales = "free_y") +
scale_x_date(date_breaks = "1 weeks",date_labels = "%d-%b",expand = c(0,0)) +
geom_vline(xintercept = as_date("2020-03-15"),size=1.5,linetype="dashed") +
labs(title = "Electricity Consumption - Distribution Companies (Consumo Eléctrico - Empresas Distribuidoras)",subtitle = "Electricity consumption measured at stations owned by distribution companies (empresas distribuidoras).",
y = element_blank(),x=element_blank(),fill="Department (Departamento)", caption = "Source: COES. Author: Luis José Zapata Bobadilla.") +
theme_gray()+
theme(axis.text.x.bottom = element_text(angle = 90,hjust=0.7)) +
ggforce::geom_mark_ellipse(aes(filter=(gwh > 4)&(DEPARTAMENTO=="LAMBAYEQUE")),fill="red")

Key Observations:

General Trends Across Departments:
- All departments show a clear decrease in electricity consumption after March 16, 2020, reflecting the impact of quarantine measures on economic activity.
- Departments like Lima and Piura, which contribute heavily to overall electricity demand, show a substantial drop.
Notable Patterns:
- The periodic spikes suggest weekly consumption cycles, particularly visible in regions like Lima and Piura.
Outlier in Lambayeque:
- The Lambayeque panel highlights an outlier (shaded red) after March 30, indicating a data registration error where recorded consumption was unrealistically high. This emphasizes the need for thorough data cleaning and verification to ensure accurate analysis.
Consistency of Patterns:
- While the drop after quarantine is expected, the consistent patterns across departments highlight the robustness of the recorded trends, aside from the noted outlier in Lambayeque, when we aggregate by points of measurement.

Addressing redudant information

A large portion of the information measured in the system is redundant, meaning it is repeated. As more measurement points (puntos de medición) are added, the redundant information in each department (departamento) increases exponentially. If demand is summed by department without applying any filters, some departments may display consumption higher than the actual values. This is the case for Lima in our dataset.

An effective way to address this problem is to exclude every measurement point (punto de medición) and instead focus only on those points that sum the majority of the measurement points in each department (departamento), as seen with the Zorritos - Tumbes point in the earlier example.

## Joining with `by = join_by(PUNTO.DE.MEDICIÓN)`
## `summarise()` has grouped output by 'Fecha'. You can override using the
## `.groups` argument.

The plot provides a visual representation of the changes in demand over time, categorized by department.

Key Observations:

Temporal Trends:
- The x-axis represents dates, ranging from mid-February to mid-April, showing daily electricity consumption.
- There is a noticeable decrease in electricity demand starting around March 16, 2020, marked by the vertical dashed line. This corresponds to the start of the COVID-19 quarantine period in Peru, reflecting reduced industrial and business activity.
Departmental Breakdown:
- The stacked area chart highlights the contribution of each department to the total electricity demand.
- Lima (pink) consistently accounts for the largest share of electricity consumption, reflecting its role as the country’s economic hub.
- Other significant contributors include Piura (blue), La Libertad (green), and Arequipa (teal).
Impact of Quarantine:
- The reduction in demand after March 16 is visible across all departments, with a sharp decline in Lima and noticeable decreases in other regions like Arequipa and Lambayeque.
- Smaller departments grouped under “OTROS” (Others) also show a slight decline, though their contribution to overall demand is relatively small.
Seasonal or Weekly Patterns:
- The wavelike structure in the plot suggests periodic fluctuations, likely driven by weekday-weekend patterns in electricity usage.

Insights:

This plot effectively demonstrates how a major event (quarantine) can impact electricity demand across different regions. It highlights the significance of Lima and other key regions while showcasing how such analyses in real time can provide actionable insights for policymakers and energy planners.

Conclusion

As mentioned earlier, the best way to avoid these types of problems is to focus solely on measurement points (puntos de medición) that aggregate the majority of other points (such as the Zorritos - Tumbes point).

To achieve this, it is necessary to have regional knowledge of which points aggregate others. This can be done by reviewing the single-line diagram (diagrama unifilar) from COES.

It is essential to identify the relevant measurement points (puntos de medición) that prevent redundancies in each region.

This is a pending task that would enable a reliable analysis by department (departamento).

It is worth noting that COES already publishes compiled information for three zones of Peru (North, Central, and South) in its Daily Evaluation Report (IEOD) under the sheet “Demand by Zone” (Demanda por zona). Therefore, if the goal is to obtain compiled information for these zones, it is easier to download it directly from the COES website.

However, if the objective is to analyze the information by department, the approach outlined in this analysis would be the most appropriate.

This analysis illustrates my proficiency in data cleaning, statistical modeling, and visualization using R, along with my ability to adapt to technical topics like electricity distribution systems.

All the files used in this analysis can be downloaded from the following link:
https://github.com/luis-zapata-b/Datos/archive/master.zip

Final Conclusion

This analysis highlights the critical role of data-driven methodologies in understanding electricity demand at the regional level, emphasizing its importance for policy-making, investment decisions, and sustainability efforts. By addressing key challenges such as redundant data, missing values, and outliers, the approach presented ensures accurate and actionable insights into energy consumption trends across Peru’s departments.

The results demonstrate how external factors, like the COVID-19 quarantine, can significantly impact electricity demand, with clear patterns observed across major regions like Lima, Piura, and Lambayeque. Moreover, the use of advanced imputation methods and outlier filtering techniques underscores the importance of rigorous data preparation in ensuring reliable conclusions.

The mechanism developed in this study not only provides valuable insights into regional electricity consumption but also establishes a reproducible framework. By leveraging publicly available data from COES, this methodology can be replicated and scaled for future analyses, supporting government efforts to optimize resource allocation and plan for sustainable energy systems.

Overall, this project showcases the integration of technical expertise in R, critical thinking, and adaptability to complex topics, demonstrating a robust foundation for tackling real-world challenges in resource management and energy analysis.

Exploring Neural Networks with Pytorch

Sat, 23 Nov 2024 00:00:00 +0000

In this blog, we’ll explore how to use a neural network built with PyTorch to predict handwritten digits from the MNIST dataset. Along the way, we’ll break down key concepts like tensors, neural network architecture, and training processes, making it easier to understand the fundamental building blocks of deep learning.

The MNIST dataset, a collection of 28x28 grayscale images of digits (0–9), is a classic benchmark for machine learning. Our goal is to create a model that accurately recognizes these digits using the power of neural networks.

Introduction

Machine learning is a transformative tool for modeling complex relationships in data. Neural networks, a class of machine learning models inspired by the human brain, are particularly adept at handling non-linear relationships. This post explores their application, focusing on tensors, neural networks, and their underlying theory.

Tensors: The Backbone of Neural Networks

At the heart of neural networks lie tensors, multidimensional arrays that generalize matrices to higher dimensions. They efficiently handle complex data like images, text, or time series, making them critical for machine learning.

For instance, a scalar is a rank-0 tensor (a single number), a vector is rank-1 (a sequence of numbers), and a matrix is rank-2 (a 2D grid). Consider This 2D input tensor x, which could represent a dataset:

import torch
x = torch.tensor([[1, 2], [3, 4], [5, 6]]) # A 2D tensor (matrix)
print(x)

## tensor([[1, 2],
## [3, 4],
## [5, 6]])

Tensors flow through the network, undergoing operations like matrix multiplication and activation functions, ultimately producing outputs. Their ability to efficiently process data on GPUs makes tensors indispensable for modern machine learning.

Neural Network Architecture

A neural network is composed of layers of interconnected nodes (neurons) that process input data and produce predictions:

Input Layer: Receives the raw data (e.g., tensors representing images, text, or numbers) and passes it to the next layer.
Hidden Layers: These intermediate layers perform computations, applying weights, biases, and activation functions to extract meaningful patterns and features.
Output Layer: Generates the final output, such as classifications, predictions, or probabilities.

Each connection between neurons has a weight that adjusts during training, enabling the network to learn from data by minimizing error. Activation functions like ReLU or Sigmoid introduce non-linearity, allowing the network to model complex relationships.

Example Model Configuration

The neural network model in this example uses a simple feedforward architecture. We define it as:

import torch.nn as nn
class SimpleNN(nn.Module):
def __init__(self, input_size, hidden_size, output_size):
super(SimpleNN, self).__init__()
self.fc1 = nn.Linear(input_size, hidden_size)
self.relu = nn.ReLU()
self.fc2 = nn.Linear(hidden_size, output_size)
def forward(self, x):
x = self.fc1(x)
x = self.relu(x)
x = self.fc2(x)
return x
input_size = 28 * 28 # Example: 28x28 images
hidden_size = 128
output_size = 10 # Example: 10 classes
model = SimpleNN(input_size, hidden_size, output_size)

Lets dive into a more detailed code explanation of the neural network architecture.

1. The Dataset: MNIST

The MNIST dataset is a collection of grayscale images of handwritten digits, commonly used as a benchmark dataset in machine learning. Each image is:

28x28 pixels in size.
Labeled with a number from 0 to 9 (10 classes in total).

Code Explanation

from torchvision import datasets, transforms
# Define transformations for the dataset
transform = transforms.Compose([transforms.ToTensor()])

Transforms: The dataset is preprocessed using transforms.ToTensor(), which converts the images into PyTorch tensors (normalized to [0, 1]).
Compose: Allows multiple transformations to be chained together. Here, only the conversion to a tensor is applied.

# Load training and testing datasets
train_dataset = datasets.MNIST(root='./data', train=True, transform=transform, download=True)
test_dataset = datasets.MNIST(root='./data', train=False, transform=transform, download=True)

train=True/False: Indicates whether to load the training or testing split of the MNIST dataset.

# Create DataLoaders
train_loader = DataLoader(train_dataset, batch_size=100, shuffle=True)
test_loader = DataLoader(test_dataset, batch_size=100, shuffle=False)

DataLoader: A utility that:
- Loads data in batches (size batch_size=100 here).
- Shuffles the training data for randomness.
- Ensures efficient data loading during training/testing.

A batch refers to a subset of the dataset processed at one time during training or inference. Instead of feeding the entire dataset to the model, which can be computationally expensive and memory-intensive, the data is divided into smaller groups (batches). Each batch passes through the network, computes predictions, and calculates the loss. The model then updates its weights based on the average loss of the batch. This approach balances efficiency and learning, making training manageable even for large datasets like MNIST.

2. The Neural Network: `SimpleNN`

A neural network is composed of layers of interconnected nodes (neurons) that process input data and produce predictions. Here’s how SimpleNN is structured:

Code Explanation

class SimpleNN(nn.Module):
def __init__(self, input_size, hidden_size, output_size):
super(SimpleNN, self).__init__()
self.fc1 = nn.Linear(input_size, hidden_size)
self.relu = nn.ReLU()
self.fc2 = nn.Linear(hidden_size, output_size)

nn.Module: Base class for all PyTorch models.
__init__: Defines the structure of the neural network.
- input_size: The number of input features (28x28 = 784 for MNIST images flattened into vectors).
- hidden_size: Number of neurons in the hidden layer (set to 128 here).
- output_size: The number of output classes (10 digits: 0 to 9).
- fc1: A fully connected (dense) layer that transforms the input into the hidden layer with dimensions [input_size -> hidden_size].
- ReLU: Applies the ReLU activation function, introducing non-linearity to the model. ReLU replaces negative values with 0.
- fc2: A fully connected layer that transforms the hidden layer output to the final output layer with dimensions [hidden_size -> output_size].

 def forward(self, x):
x = self.fc1(x)
x = self.relu(x)
x = self.fc2(x)
return x

forward: Defines how data flows through the network.

Passes the input through fc1 (input layer).
Applies the ReLU activation.
Passes the result through fc2 (output layer).
Returns the final output.

init establishes the components of the network, and forward dictates how those components are used during the forward pass. Together, they form the blueprint and operational flow of the model.

The ReLU activation function (Rectified Linear Unit) is one of the most commonly used activation functions in neural networks, especially in deep learning. Its popularity stems from its simplicity, efficiency, and effectiveness.

What Is ReLU?

The ReLU function is defined mathematically as:

$$ f(x) = \begin{cases} x & \text{if } x > 0 \\
0 & \text{if } x \leq 0 \end{cases} $$ Or simply:

$$ f(x) = \max(0, x) $$ In words, ReLU outputs the input value directly if it’s positive and outputs 0 for any negative value.

Why Use an Activation Function?

Activation functions introduce non-linearity into the neural network. Without an activation function, the network would simply consist of linear transformations (matrix multiplications), which means it would only be able to model linear relationships in data. Non-linear activation functions, like ReLU, allow the network to model complex, non-linear relationships.

ReLU is applied to the output of each neuron (or unit) in a layer after the weighted sum of inputs and bias. Here’s what happens:

A neuron calculates a weighted sum of its inputs:

$$ z = \mathbf{w} \cdot \mathbf{x} + b, \\ \
where: \\ \
\mathbf{w}: Weights, \\ \
\mathbf{x}: Inputs, \\ \
b: Bias $$

The ReLU function is applied to (z) to produce the neuron’s output:

$$a = f(z) = \max(0, z)$$

This process is repeated across all neurons in the layer.

Benefits of ReLU

Simplicity: ReLU is computationally inexpensive because it involves simple thresholding at zero. Unlike other activation functions, like sigmoid or tanh, it does not require complex mathematical operations like exponentials.
Avoids Vanishing Gradients:
- Activation functions like sigmoid and tanh compress outputs into a small range, which can cause gradients to vanish during backpropagation, making learning very slow in deep networks.
- ReLU doesn’t compress values, so the gradient is preserved for positive inputs, enabling faster and more stable learning.
Efficient Sparsity:
- Since ReLU outputs 0 for negative inputs, it introduces sparsity in the network. Sparsity helps by reducing the dependency between neurons and improving computation efficiency.

ReLU introduces non-linearity between the input and output layers, allowing the network to learn complex patterns in the MNIST dataset.
Its simplicity and efficiency make it an ideal choice for a small-scale network like this one.

Visualizing ReLU

Here’s how ReLU behaves visually:

For $(x > 0)$, ReLU outputs the same value.
For $(x \leq 0)$, ReLU outputs 0.

In summary, ReLU is simple, fast, and effective for training deep neural networks, making it the default choice in many modern architectures.

3. Model Initialization

input_size = 28 * 28 # Example: 28x28 images
hidden_size = 128
output_size = 10 # Example: 10 classes
model = SimpleNN(input_size, hidden_size, output_size)

input_size = 28 * 28: Each MNIST image is a 28x28 matrix. Before feeding it into the model, it’s flattened into a 1D tensor of size 784.
hidden_size = 128: The hidden layer contains 128 neurons. This is a hyperparameter and can be adjusted.
output_size = 10: The output layer contains 10 neurons, one for each digit class (0 to 9).
model: Instantiates the SimpleNN class with the specified parameters.

4. Model Training Expectations

The goal is to train the model to correctly classify MNIST digits with high accuracy. Here’s what we expect:

Input: A batch of flattened MNIST images ([batch_size, 784]).
Hidden Layer: Maps the input to a feature space of size 128.
Output Layer: Produces probabilities for each digit class using a softmax layer (applied externally during training).
Loss Function: Measures how far the predictions are from the true labels.
Optimization: Uses Stochastic Gradient Descent (SGD) or another optimizer to minimize the loss by adjusting weights.

The loss function serves as a key metric in machine learning, quantifying how well or poorly a model’s predictions align with the actual target values. By evaluating the difference between the predicted outcomes and the true labels, the loss function provides a numerical value that guides the optimization process, helping to improve the model’s performance.

To minimize this loss, gradient descent is often used as an optimization algorithm. It works by iteratively adjusting the model’s parameters in the direction that reduces the loss function the most. By following the gradient (or slope) of the loss function, gradient descent ensures the model progressively improves its predictions, moving closer to the optimal solution.

5. Expected Results

Accuracy: For a simple feedforward network like this, we expect around 90%-92% accuracy on the MNIST test set after training.
Loss: The loss will start high (e.g., 2.0) and gradually decrease to values around 0.3-0.5 as the model learns.
Improvement Potential:
- Adding more hidden layers or increasing the hidden size.
- Experimenting with advanced architectures like CNNs.
- Using data augmentation techniques for better generalization.

Example

Lets combine this concepts, and run the model in order the get the results.

Data Loading

We use PyTorch’s DataLoader to load the training and testing data.

from torch.utils.data import DataLoader
from torchvision import datasets, transforms
# Define transformations for the dataset
transform = transforms.Compose([transforms.ToTensor()])
# Load training and testing datasets
train_dataset = datasets.MNIST(root='./data', train=True, transform=transform, download=True)
test_dataset = datasets.MNIST(root='./data', train=False, transform=transform, download=True)
# Create DataLoaders
train_loader = DataLoader(train_dataset, batch_size=100, shuffle=True)
test_loader = DataLoader(test_dataset, batch_size=100, shuffle=False)

Loss Function and Optimizer

We use Cross Entropy Loss as the loss function and SGD as the optimizer:

import torch.optim as optim
loss_fn = nn.CrossEntropyLoss()
optimizer = optim.SGD(model.parameters(), lr=0.01)

Training the Neural Network

Training a neural network involves repeatedly feeding data through the model, calculating the error (loss), and updating the model’s parameters using optimization techniques like gradient descent. The train_loop() function handles this process for each training epoch.

Here’s how train_loop() is implemented:

def train_loop(dataloader, model, loss_fn, optimizer):
"""
Training loop for a neural network model.
Parameters:
- dataloader: DataLoader object for training data
- model: Neural network model
- loss_fn: Loss function (e.g., CrossEntropyLoss)
- optimizer: Optimizer (e.g., SGD, Adam)
Returns:
None
"""
size = len(dataloader.dataset)
model.train()
for batch, (X, y) in enumerate(dataloader):
# Flatten images into vectors
X = X.view(X.size(0), -1)
# Forward pass
pred = model(X)
loss = loss_fn(pred, y)
# Backpropagation
optimizer.zero_grad()
loss.backward()
optimizer.step()
# Log progress
if batch % 100 == 0:
loss_current = loss.item()
print(f"loss: {loss_current:>7f} [{batch * len(X):>5d}/{size:>5d}]")

This function is executed within each training epoch, as shown below:

for epoch in range(5):
train_loop(train_loader, model, loss_fn, optimizer)

## loss: 2.302421 [ 0/60000]
## loss: 2.082567 [10000/60000]
## loss: 1.820236 [20000/60000]
## loss: 1.481050 [30000/60000]
## loss: 1.281981 [40000/60000]
## loss: 1.009585 [50000/60000]
## loss: 0.945703 [ 0/60000]
## loss: 0.837330 [10000/60000]
## loss: 0.799911 [20000/60000]
## loss: 0.746439 [30000/60000]
## loss: 0.629169 [40000/60000]
## loss: 0.529109 [50000/60000]
## loss: 0.390061 [ 0/60000]
## loss: 0.420087 [10000/60000]
## loss: 0.353291 [20000/60000]
## loss: 0.467722 [30000/60000]
## loss: 0.710849 [40000/60000]
## loss: 0.349696 [50000/60000]
## loss: 0.385319 [ 0/60000]
## loss: 0.443223 [10000/60000]
## loss: 0.247738 [20000/60000]
## loss: 0.438087 [30000/60000]
## loss: 0.547100 [40000/60000]
## loss: 0.468050 [50000/60000]
## loss: 0.330704 [ 0/60000]
## loss: 0.365649 [10000/60000]
## loss: 0.326305 [20000/60000]
## loss: 0.363884 [30000/60000]
## loss: 0.315488 [40000/60000]
## loss: 0.360471 [50000/60000]

Visualizing Predictions

After training, we can compare the model’s predictions to true labels using a histogram:

import matplotlib.pyplot as plt
def test_loop(dataloader, model, loss_fn):
model.eval()
size = len(dataloader.dataset)
num_batches = len(dataloader)
test_loss = 0
correct = 0
all_preds = []
all_labels = []
with torch.no_grad():
for X, y in dataloader:
X = X.view(X.size(0), -1)
pred = model(X)
test_loss += loss_fn(pred, y).item()
correct += (pred.argmax(1) == y).type(torch.float).sum().item()
all_preds.extend(pred.argmax(1).tolist())
all_labels.extend(y.tolist())
test_loss /= num_batches
correct /= size
return test_loss, correct, all_preds, all_labels
test_loss, test_acc, pred, lab = test_loop(test_loader, model, loss_fn)
print(f"Test Error: \n Accuracy: {(100*test_acc):>0.1f}%, Avg loss: {test_loss:>8f} \n")

## Test Error:
## Accuracy: 90.5%, Avg loss: 0.341179

# Plot histograms
plt.hist(pred, bins=10, alpha=0.5, label='Predictions')
plt.hist(lab, bins=10, alpha=0.5, label='Labels')
plt.legend()
plt.show()

Conclusion

In this project, we trained a simple feedforward neural network on the MNIST dataset and evaluated its performance. The final results were as follows:

Test Accuracy: 90.4%
Average Test Loss: 0.3445

The histogram above shows that the predictions closely match the true labels, demonstrating that the neural network successfully learned patterns in the data. While the accuracy indicates solid performance, further improvements could be achieved through:

Optimizing hyperparameters (e.g., learning rate, batch size, number of hidden units).
Experimenting with deeper or more complex architectures like Convolutional Neural Networks (CNNs).
Increasing the number of training epochs to allow the network to converge better.

This blog demonstrated how to build, train, and evaluate a simple feedforward neural network using PyTorch to classify handwritten digits from the MNIST dataset. Along the way, we covered essential concepts such as tensors, neural network architecture, activation functions, and the training process.

Neural networks are versatile tools with applications across diverse domains. This project serves as a foundation for understanding how to implement and optimize these models, paving the way for tackling more complex machine learning challenges in the future.

References

PyTorch Documentation: https://pytorch.org/docs/
MNIST Dataset: http://yann.lecun.com/exdb/mnist/

Kernel Density and Local Estimators with Confidence Bands

Sat, 23 Nov 2024 00:00:00 +0000

Introduction

In this post, we’ll explore Kernel Density Estimation (KDE), Local Constant Estimator, and Local Linear Estimator. Additionally, we’ll introduce 95% Confidence Bands (both point-wise and universal) to quantify the uncertainty in the estimators.

This will provide a more complete understanding of these non-parametric methods and how to use them effectively.

Kernel Density Estimation (KDE) with Confidence Bands

Kernel Density Estimation (KDE) estimates the probability density function of a random variable by smoothing data points with a kernel function. To evaluate the reliability of the KDE, we add point-wise and universal confidence bands.

Kernel Density Estimation (KDE): Theory and Application

Kernel Density Estimation (KDE) is a non-parametric method used to estimate the probability density function (PDF) of a dataset. Unlike parametric approaches, KDE does not assume the data follows a specific distribution (e.g., normal distribution). Instead, it constructs the density estimate by smoothing individual data points using a kernel function. This makes KDE highly flexible and particularly useful when the underlying data distribution is unknown or non-standard, such as in financial data where heavy tails and skewness are common.

The KDE is defined as:

$$ \hat{f}(x) = \frac{1}{n h} \sum_{i=1}^n K\left(\frac{x - X_i}{h}\right), $$

where:

$ \hat{f}(x) $ is the estimated density at point $ x $,
$ X_i $ are the data points,
$ K $ is the kernel function (e.g., Epanechnikov or Gaussian),
$ h $ is the bandwidth, controlling the smoothness of the estimate.

The choice of bandwidth $ h $ is crucial. Smaller $ h $ values lead to a more flexible estimate but can result in overfitting, while larger $ h $ values produce a smoother density at the cost of losing detail.

When to Use KDE

KDE is particularly valuable in situations where the density of the data is unknown and may deviate from standard distributions. For example:

Exploratory Analysis: When we want to visualize the shape of the data’s distribution without assuming normality.
Non-Normal Data: Financial data, for instance, often exhibit heavy tails or skewness, making KDE a robust choice to uncover the true density.

In the code example above, the KDE is compared to a normal density estimation.

Confidence Bands

Confidence bands, both pointwise and uniform, provide uncertainty estimates for the density function:

Pointwise Confidence Bands quantify the variability at specific points, showing how confident we are about the estimate at each ( x ).
Uniform Confidence Bands cover the entire range of ( x ), ensuring that the true density lies within the band with a given probability.

These bands help evaluate the reliability of the KDE and highlight areas where the estimate may be less accurate due to data sparsity or noise.

Example: KDE with Confidence Bands

Here’s a KDE estimator example with pointwise confidence bands for financial data using the Epanechnikov kernel. The code downloads stock price data (APPLE), computes log returns, and estimates the density with confidence bands. The normal density is also calculated for comparison. :

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import yfinance as yf
# Download financial data
def download_financial_data(ticker, start_date, end_date):
data = yf.download(ticker, start=start_date, end=end_date)
log_returns = np.log(data['Close'] / data['Close'].shift(1)).dropna()
return log_returns
# Define KDE with Epanechnikov Kernel
def epanechnikov_kernel(u):
return np.where(np.abs(u) <= np.sqrt(5), (3 / (4 * np.sqrt(5))) * (1 - u**2 / 5), 0)
def kernel_density_estimate(x_values, data, bandwidth):
densities = np.zeros_like(x_values)
for i, x in enumerate(x_values):
weights = epanechnikov_kernel((x - data) / bandwidth)
densities[i] = np.sum(weights) / (len(data) * bandwidth)
return densities
# Function to calculate variance estimate for pointwise confidence intervals
def variance_estimate(x_values, data, bandwidth):
n = len(data)
variances = np.zeros_like(x_values)
for i, x0 in enumerate(x_values):
K_h_values = epanechnikov_kernel((data - x0) / bandwidth) / bandwidth
term1 = (1 / n) * np.sum(K_h_values**2)
term2 = ((1 / n) * np.sum(K_h_values))**2
variances[i] = (term1 - term2) / n
return variances
# Function to calculate the normal density
def normal_density_estimate(x_values, mean, std_dev):
return (1 / (np.sqrt(2 * np.pi * std_dev**2))) * np.exp(-((x_values - mean)**2) / (2 * std_dev**2))
# Download and preprocess financial data
ticker = "AAPL" # Example ticker
start_date = "2023-01-01"
end_date = "2024-01-01"
log_returns = download_financial_data(ticker, start_date, end_date)

## /Users/luisj/anaconda3/lib/python3.11/site-packages/yfinance/utils.py:775: FutureWarning: The 'unit' keyword in TimedeltaIndex construction is deprecated and will be removed in a future version. Use pd.to_timedelta instead.
## df.index += _pd.TimedeltaIndex(dst_error_hours, 'h')
##
[*********************100%%**********************] 1 of 1 completed

# KDE and confidence bands
x_values = np.linspace(log_returns.min(), log_returns.max(), 1000)
bandwidth = 0.02 # Example bandwidth
density_estimates = kernel_density_estimate(x_values, log_returns.values, bandwidth)
# Pointwise confidence bands
alpha = 0.05
variance_estimates = variance_estimate(x_values, log_returns.values, bandwidth)
z_alpha = norm.ppf(1 - alpha / 2)
lower_bounds = density_estimates - z_alpha * np.sqrt(variance_estimates)
upper_bounds = density_estimates + z_alpha * np.sqrt(variance_estimates)
# Normal density for comparison
mean_returns = log_returns.mean()
std_returns = log_returns.std()
normal_density = normal_density_estimate(x_values, mean_returns, std_returns)
# Plot results
plt.figure(figsize=(12, 6))
plt.plot(x_values, density_estimates, label="KDE (Epanechnikov)", color='blue')
plt.plot(x_values, normal_density, label="Normal Density (Mean & Std Dev)", color='red', linestyle='--')
plt.fill_between(x_values, lower_bounds, upper_bounds, color='gray', alpha=0.3, label="95% Pointwise Confidence Band")
plt.xlabel("Log Returns")
plt.ylabel("Density")
plt.title(f"Kernel Density Estimation for {ticker} with Confidence Bands")
plt.legend()
plt.show()

Explanation:

Epanechnikov Kernel: The function smooths the data points based on the Epanechnikov formula. A bandwidth parameter controls the degree of smoothing.
Normal Density: The normal density estimate is calculated using the sample mean and standard deviation for comparison with the KDE.
Confidence Bands:
- Variance Estimate: Variance for pointwise confidence intervals is computed based on kernel weights.

Conclusion

KDE is a powerful tool for density estimation, especially when the underlying distribution of the data is unknown or non-standard. It is particularly useful in domains like finance, where data often deviate from normality, exhibiting heavy tails or irregular shapes. By using kernels and bandwidths, KDE provides a flexible, data-driven approach to estimate densities without relying on strict parametric assumptions. The addition of confidence bands further enhances its utility by quantifying the uncertainty in the estimates, making KDE a robust choice for both exploratory and confirmatory data analysis.

Introduction: Nonparametric Regression

Using the estimation of Kernel Density Estimators, now we can approach the expected values of Y given X (using the esitmated density). By doing this, we can obtain non parametric estimators. Nonparametric regression methods are powerful tools for estimating the relationship between a response variable $ Y $ and a predictor $ X $, without assuming a specific parametric form for the underlying relationship. Unlike traditional linear or polynomial regression, nonparametric methods adapt to the data’s structure, providing flexibility in modeling complex relationships.

Two widely used nonparametric regression estimators are the Nadaraya-Watson Estimator and the Local Linear Estimator. While both methods use kernel functions to assign weights to data points based on their proximity to the evaluation point, they differ in how they model the local behavior of the data.

The Nadaraya-Watson Estimator

The Nadaraya-Watson Estimator is a simple, kernel-based approach to regression. For a given point $x $, the estimator computes the weighted average of the response variable $ Y $ using the weights determined by the kernel function:

$$ \hat{m}(x) = \frac{\sum_{i=1}^n K_h(x - X_i) Y_i}{\sum_{i=1}^n K_h(x - X_i)}, $$

where:

$( K_h(u) = \frac{1}{h} K\left(\frac{u}{h}\right) ) $ is the kernel function scaled by the bandwidth ( h ),
$( X_i )$ and $( Y_i )$ are the predictor and response variables, respectively.

The kernel function ( K ) determines the weight assigned to each observation. Common kernels include the Gaussian kernel and the Epanechnikov kernel. The bandwidth ( h ) controls the smoothness of the estimate: smaller values of ( h ) result in a more flexible, less smooth fit, while larger values yield smoother estimates.

While the Nadaraya-Watson Estimator is simple and effective, it can suffer from boundary bias. Near the edges of the data, the kernel weights become asymmetric, leading to inaccuracies in the estimated values.

The Local Linear Estimator

The Local Linear Estimator extends the Nadaraya-Watson approach by fitting a local linear regression in the neighborhood of each evaluation point ( x ). Instead of assuming a constant value for ( Y ) (as in NW), this method models ( Y ) as a linear function:

$$ Y \approx \theta_0 + \theta_1 (X - x).$$

For a given ( x ), the parameters $ \theta_0 $ and $ \theta_1 $ are estimated by minimizing a kernel-weighted least-squares objective:

$$ \min_{\theta_0, \theta_1} \sum_{i=1}^n K_h(X_i - x) \left[ Y_i - \theta_0 - \theta_1 (X_i - x) \right]^2. $$

The solution gives the local linear regression coefficients $ \hat{\theta}_0 $ and $ \hat{\theta}_1 $, and the estimator at ( x ) is $ \hat{m}(x) = \hat{\theta}_0 $.

By allowing for a linear fit within the kernel window, the Local Linear Estimator addresses the boundary bias issue present in the Nadaraya-Watson Estimator. It can better capture the local structure of the data and typically performs well both near the edges and in regions with varying density.

Comparison of the Two Methods

Feature	Nadaraya-Watson Estimator	Local Linear Estimator
Model Assumption	Locally constant	Locally linear
Bias near boundaries	High (prone to boundary bias)	Low (reduces boundary bias)
Complexity	Simple (weighted average)	More complex (weighted least squares)
Behavior near edges	Underestimates or overestimates response	Handles boundaries more accurately

Both methods depend heavily on the choice of the kernel function and the bandwidth parameter, which control the smoothness and flexibility of the estimates.

Python Code

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.integrate import quad
# Define kernels
def epanechnikov_kernel(u):
return np.where(np.abs(u) <= np.sqrt(5), (3 / (4 * np.sqrt(5))) * (1 - u**2 / 5), 0)
# Nadaraya-Watson Estimator
def nadaraya_watson(X, Y, x_grid, bandwidth, kernel=epanechnikov_kernel):
estimates = []
for x in x_grid:
weights = kernel((X - x) / bandwidth)
numerator = np.sum(weights * Y)
denominator = np.sum(weights)
estimates.append(numerator / denominator if denominator > 0 else 0)
return np.array(estimates)
# Local Linear Estimator
def local_linear_estimator(X, Y, x_eval, bandwidth, kernel=epanechnikov_kernel):
n = len(X)
U = np.vstack([np.ones(n), X - x_eval]).T
W = np.diag(kernel((X - x_eval) / bandwidth) / bandwidth)
theta_hat = np.linalg.inv(U.T @ W @ U) @ (U.T @ W @ Y)
return theta_hat[0], theta_hat[1]
# Variance estimator for pointwise confidence intervals
def compute_s_hat(X, residuals, x_eval, bandwidth, kernel=epanechnikov_kernel):
weights = kernel((X - x_eval) / bandwidth) / bandwidth
weights_squared = weights**2
residuals_squared = residuals**2
f_hat = np.sum(weights) / (len(X) * bandwidth)
s_hat = np.sqrt(np.sum(weights_squared * residuals_squared) / (len(X)**2 * f_hat**2))
return s_hat
# Generate synthetic data
np.random.seed(42)
X = np.linspace(0, 10, 100)
Y = np.sin(X) + np.random.normal(0, 0.2, size=len(X))
# Define evaluation grid and bandwidth
x_grid = np.linspace(0, 10, 200)
bandwidth = 0.8
# Estimate Nadaraya-Watson
nw_estimates = nadaraya_watson(X, Y, x_grid, bandwidth)
# Estimate Local Linear
ll_intercepts = []
for x in x_grid:
intercept, _ = local_linear_estimator(X, Y, x, bandwidth)
ll_intercepts.append(intercept)
# Residuals for Nadaraya-Watson
nw_residuals = Y - nadaraya_watson(X, Y, X, bandwidth)
# Pointwise Confidence Intervals
alpha = 0.05
z_alpha = norm.ppf(1 - alpha / 2)
pointwise_lower = []
pointwise_upper = []
s_hat_values = []
for x_eval, m_eval in zip(x_grid, nw_estimates):
s_hat = compute_s_hat(X, nw_residuals, x_eval, bandwidth)
s_hat_values.append(s_hat)
pointwise_lower.append(m_eval - z_alpha * s_hat)
pointwise_upper.append(m_eval + z_alpha * s_hat)
# Plot Nadaraya-Watson with Pointwise CI
plt.figure(figsize=(12, 6))

## <Figure size 1200x600 with 0 Axes>

plt.plot(x_grid, nw_estimates, label="Nadaraya-Watson Estimate", color="blue")

## [<matplotlib.lines.Line2D object at 0x17d060550>]

plt.fill_between(x_grid, pointwise_lower, pointwise_upper, color="blue", alpha=0.2, label="95% Pointwise CI")

## <matplotlib.collections.PolyCollection object at 0x17d0f1b90>

plt.scatter(X, Y, alpha=0.5, label="Data Points", color="black", s=10)

## <matplotlib.collections.PathCollection object at 0x17d0e6cd0>

plt.title("Nadaraya-Watson Estimator with Pointwise Confidence Intervals")

## Text(0.5, 1.0, 'Nadaraya-Watson Estimator with Pointwise Confidence Intervals')

plt.xlabel("X")

## Text(0.5, 0, 'X')

plt.ylabel("Y")

## Text(0, 0.5, 'Y')

plt.legend()

## <matplotlib.legend.Legend object at 0x17d069d10>

plt.show()

# Uniform Confidence Bands
def compute_G(x_grid, X, residuals, bandwidth, kernel=epanechnikov_kernel):
n = len(X)
G_values = []
for x_eval in x_grid:
kernel_values = kernel((X - x_eval) / bandwidth) / bandwidth
total_weight = np.sum(kernel_values)
weights = kernel_values / total_weight if total_weight > 0 else np.zeros_like(kernel_values)
xi = np.random.normal(0, 1, n)
s_hat = compute_s_hat(X, residuals, x_eval, bandwidth)
G_x = np.sum(xi * residuals * weights / s_hat)
G_values.append(G_x)
return np.array(G_values)
def gaussian_multiplier_bootstrap(x_grid, X, residuals, bandwidth, kernel=epanechnikov_kernel, B=1000, alpha=0.05):
bootstrap_sup = []
for _ in range(B):
G_values = compute_G(x_grid, X, residuals, bandwidth, kernel)
bootstrap_sup.append(np.max(np.abs(G_values)))
return np.quantile(bootstrap_sup, 1 - alpha)
q_alpha = gaussian_multiplier_bootstrap(x_grid, X, nw_residuals, bandwidth)
uniform_lower = nw_estimates - q_alpha * np.array(s_hat_values)
uniform_upper = nw_estimates + q_alpha * np.array(s_hat_values)
# Plot with Uniform Confidence Bands
plt.figure(figsize=(12, 6))

## <Figure size 1200x600 with 0 Axes>

plt.plot(x_grid, nw_estimates, label="Nadaraya-Watson Estimate", color="blue")

## [<matplotlib.lines.Line2D object at 0x17de5ebd0>]

plt.fill_between(x_grid, uniform_lower, uniform_upper, color="red", alpha=0.1, label="95% Uniform CI")

## <matplotlib.collections.PolyCollection object at 0x17de92590>

plt.scatter(X, Y, alpha=0.5, label="Data Points", color="black", s=10)

## <matplotlib.collections.PathCollection object at 0x17de2c410>

plt.title("Nadaraya-Watson Estimator with Uniform Confidence Bands")

## Text(0.5, 1.0, 'Nadaraya-Watson Estimator with Uniform Confidence Bands')

plt.xlabel("X")

## Text(0.5, 0, 'X')

plt.ylabel("Y")

## Text(0, 0.5, 'Y')

plt.legend()

## <matplotlib.legend.Legend object at 0x17de9ee90>

plt.show()

Pointwise vs. Uniform Confidence Bands

Pointwise and uniform confidence bands both quantify the uncertainty in regression estimates but differ in their scope. Pointwise confidence bands provide an interval estimate for the regression function at each specific point ( x ), indicating the range where the true function is likely to lie with a given probability (e.g., 95%). However, these intervals do not account for the simultaneous uncertainty across all points in the domain.

In contrast, uniform confidence bands account for the entire range of ( x ), ensuring that the true regression function lies within the band at all points with a specified probability (e.g., 95%). While pointwise bands are narrower and easier to compute, uniform bands are more conservative and robust, capturing the overall uncertainty across the entire function.

Key Features

Synthetic Data: Non-linear function $ Y = \sin(X) + \text{noise} $.
Nadaraya-Watson Estimation: Smoothed estimates with Epanechnikov kernel.
Local Linear Estimation: Estimates the function with reduced bias.
Pointwise Confidence Intervals: Calculated using the standard error of residuals.
Uniform Confidence Bands: Computed via Gaussian multiplier bootstrap.

This approach demonstrates the application of kernel-based regression estimators and provides clear visual comparisons of pointwise and uniform confidence intervals.

# Plot Local Linear Estimator
plt.figure(figsize=(12, 6))

## <Figure size 1200x600 with 0 Axes>

plt.plot(x_grid, ll_intercepts, label="Local Linear Estimate", color="green")

## [<matplotlib.lines.Line2D object at 0x17ded4890>]

plt.scatter(X, Y, alpha=0.5, label="Data Points", color="black", s=10)

## <matplotlib.collections.PathCollection object at 0x17dec5990>

plt.title("Local Linear Estimator")

## Text(0.5, 1.0, 'Local Linear Estimator')

plt.xlabel("X")

## Text(0.5, 0, 'X')

plt.ylabel("Y")

## Text(0, 0.5, 'Y')

plt.legend()

## <matplotlib.legend.Legend object at 0x17ded64d0>

plt.show()

# Plot Comparison: Nadaraya-Watson vs Local Linear Estimators
plt.figure(figsize=(12, 6))

## <Figure size 1200x600 with 0 Axes>

plt.plot(x_grid, nw_estimates, label="Nadaraya-Watson Estimate", color="blue")

## [<matplotlib.lines.Line2D object at 0x17dec8590>]

plt.plot(x_grid, ll_intercepts, label="Local Linear Estimate", color="green", linestyle="--")

## [<matplotlib.lines.Line2D object at 0x17deabfd0>]

plt.scatter(X, Y, alpha=0.5, label="Data Points", color="black", s=10)

## <matplotlib.collections.PathCollection object at 0x17dec4290>

plt.title("Comparison: Nadaraya-Watson vs Local Linear Estimators")

## Text(0.5, 1.0, 'Comparison: Nadaraya-Watson vs Local Linear Estimators')

plt.xlabel("X")

## Text(0.5, 0, 'X')

plt.ylabel("Y")

## Text(0, 0.5, 'Y')

plt.legend()

## <matplotlib.legend.Legend object at 0x17de7f210>

plt.show()

Conclusion

The plot compares the Nadaraya-Watson Estimator (solid blue line) and the Local Linear Estimator (dashed green line) against the observed data points (gray dots). Here are some key observations:

Overall Fit:
- Both estimators capture the non-linear structure of the data effectively, reflecting the sinusoidal trend of the underlying function.
- The estimators generally align well in regions where the data density is higher.
Boundary Behavior:
- At the edges of the data (e.g., near $ X = 0 $ and $ X = 10 $), the Nadaraya-Watson Estimator shows a stronger bias, curving upward or downward.
- The Local Linear Estimator performs better at the boundaries, maintaining a smoother and more accurate extrapolation due to its local linear adjustment.
Flexibility vs. Bias:
- The Nadaraya-Watson Estimator, being locally constant, is slightly less responsive to rapid changes in the data, particularly at the boundaries where it suffers from higher bias.
- The Local Linear Estimator, with its linear adjustment, is better equipped to handle such variations, leading to reduced bias and improved accuracy near the edges.
Central Region:
- In the central region of the domain (e.g., between $ X = 3 $ and $ X = 7 $), both estimators perform similarly, as the data density is sufficient for accurate local smoothing.
Practical Implications:
- The Nadaraya-Watson Estimator may be more suitable for datasets where boundary effects are not a concern and the goal is simplicity, or we are unsure of the linear assumptions in the data. Its locally constant approach avoids making strong assumptions about the underlying structure, making it a safer choice for exploratory analysis in small datasets
- The Local Linear Estimator is preferable when boundary accuracy is critical or when the data exhibits sharp trends that require more adaptive modeling.It can handle variations better, especially when there is enough data to justify the assumption of local linearity.

Final Conclusion

In this post, we explored the theory and implementation of Kernel Density Estimation (KDE), the Nadaraya-Watson Estimator, and the Local Linear Estimator, along with methods to construct 95% confidence bands (pointwise and uniform). These non-parametric techniques provide a flexible framework for analyzing data without assuming a specific parametric model, making them powerful tools for understanding complex relationships.

Kernel Density Estimation (KDE) allows us to estimate the probability density function of a variable, with the kernel function and bandwidth playing key roles in balancing bias and variance. The inclusion of confidence bands helps quantify the uncertainty in the density estimates, providing deeper insights into the reliability of the results.

The Nadaraya-Watson Estimator is a simple and intuitive approach to nonparametric regression that uses kernel-weighted averages to estimate the response variable. While effective in capturing trends, it is susceptible to boundary bias and less responsive to sharp changes in the data.

The Local Linear Estimator, on the other hand, fits a locally linear model to reduce bias at the boundaries and adapt better to variations in the data. Its ability to handle local trends makes it a more robust choice in many scenarios, especially when data density varies across the domain.

Confidence bands enhance these methods by providing uncertainty estimates. Pointwise confidence bands quantify uncertainty at individual points, while uniform confidence bands ensure coverage over the entire domain, offering a more robust perspective on the estimator’s reliability.

The Nadaraya-Watson Estimator is particularly useful for exploratory analysis, especially with small datasets or when there are uncertainties about the linearity of the underlying relationships. In contrast, the Local Linear Estimator is better suited for situations where accuracy near boundaries or sharp data trends is critical.

In conclusion, these nonparametric methods offer flexibility and adaptability in data analysis, making them essential tools for uncovering complex patterns and relationships. By carefully selecting the appropriate estimator, kernel function, and bandwidth, you can tailor these approaches to meet your analytical objectives effectively.

Presidential Elections in Peru 2021: A Text Analysis

Fri, 22 Nov 2024 00:00:00 +0000

Introduction

How do political parties craft their messages to win over voters? This blog dives into the government plans from Peru’s 2021 presidential election, analyzing the priorities and sentiments of six major parties using cutting-edge text-mining tools. By uncovering patterns in language and tone, this analysis provides a data-driven lens into the strategies of political discourse.

This methodology could be effectively applied in other countries, like the United States, to analyze political platforms, legislative texts, or policy proposals. By extracting key terms and sentiments, it can provide insights into thematic focuses, public priorities, and ideological differences, enabling data-driven comparisons across political parties or regions.

This document analyzes the government plans for the 2021 Peruvian presidential election, plans registered with the JNE by the respective parties (the government plans can be downloaded from here). Due to the large number of parties (23), only the top 6 parties with the highest voter intentions according to IPSOS will be analyzed. The project can be replied in this Github repository: PlanGobierno.

The processing of the survey follows the methodology presented in the book Tidy Text by Julia Silge and David Robinson.

First, this analysis focuses on the government plans of the top six parties with the highest voter intentions, as reported by IPSOS: Victoria Nacional, Fuerza Popular, Partido Morado, Juntos por el Perú, Podemos Perú, and Acción Popular. Notably, Peru Libre, the party that ultimately won the election with Pedro Castillo as its candidate, was not even included in this survey. This omission starkly underscores how isolated Lima and its media are from the realities and priorities of rural Peru. Even the author, living in Lima at the time, was embedded in this bubble of perspectives, which largely ignored the growing influence of political movements outside the capital. This analysis, therefore, serves as a reminder of the biases inherent in focusing solely on urban, mainstream narratives while neglecting the broader political landscape.

The purpose of this analysis is to help the reader understand the differences between the various government plans. For this purpose, the analysis strives to present the government plans impartially.

Global Results

Below is a word cloud representing the combined text from the six government plans:

The ten most repeated words demonstrate the primary political parties’ focus on the development (desarrollo) of the country, with a strong emphasis on health (salud) in line with the current health crisis. The proposal to improve the system/services/government (sistema/servicios/gobierno) aligns with the repetition of these terms.

## # A tibble: 10 × 2
## word n
## <chr> <int>
## 1 nacional 788
## 2 desarrollo 661
## 3 salud 595
## 4 sistema 554
## 5 servicios 432
## 6 gobierno 408
## 7 social 352
## 8 plan 339
## 9 recursos 291
## 10 calidad 290

Global Sentiments

Certain words are associated with negative sentiments (e.g., problemas [problems], pobreza [poverty], crisis [crisis], etc.) and positive sentiments (e.g., reforma [reform], propuesta [proposal], bienestar [well-being], etc.). By counting these words, we can infer which sentiment is associated with the overall discourse:

Negative sentiments are observed in words such as pobreza [poverty], crisis [crisis], and lucha [struggle]. Similarly, positive sentiments are present in words like reforma [reform], derecho [rights], and capacidad [capacity].

We can infer that the government plans generally position their discourse against crisis/poverty (crisis/pobreza), offering reforms/improvements (reformas/mejoras) focused on rights/capacities (derechos/capacidades).

Comparison Between Parties

Now we will compare the government plans of each party. Below is a word cloud representing the most relevant words for each of the main parties.

Total Word Comparison

When looking at the most repeated words in each government plan, it is evident that the word salud [health] ranks among the top in all plans (except in the plan from Victoria Nacional, which emphasizes introductory terms more heavily).

The following positions provide more heterogeneous ideas about the themes relevant to each plan. It is worth mentioning that the following words were removed for being related to the general topic of governance: Nacional [National], Sistema [System], Servicios [Services], Gobierno [Government], Programa [Program], Propuesta [Proposal], Plan [Plan], and Servicios [Services]. Additionally, the following words were removed for lacking contextual relevance: Millones [Millions], Niveles [Levels], Enfoque [Focus], Forma [Form], Proponemos [We Propose], Inclusive [Inclusive], and Impulsaremos [We Will Promote].

As mentioned earlier, the most repeated words (Nacional [National], Sistema [System], Desarrollo [Development], Servicios [Services], Gobierno [Government]) are common across all parties and have been excluded. These words are all related to a shared theme within the documents: A government plan for a country. However, topics related to Salud [Health] and Desarrollo [Development] are common across all plans due to the context of COVID-19 and the economic crisis following the pandemic.

An interesting approach might be to extract the most repeated words within a single government plan that are not as commonly repeated across the other plans.

TF-IDF Comparison

The following approach uses tf-idf, a widely-used technique in information retrieval and text analysis to identify the most distinctive and relevant terms in a document.

TF (Term Frequency): Measures how often a word appears in a specific document. Words that occur frequently in a document have higher TF values.
IDF (Inverse Document Frequency): Weighs how unique a word is across an entire collection of documents. Words that appear in many documents (e.g., “government” in this context) have a lower IDF score, as they are less distinctive.

By combining these measures, tf-idf highlights words that are important to a single document but not common across the collection. This makes it ideal for comparing government plans, as it allows us to uncover terms that are particularly relevant to each party’s proposals while filtering out generic or overly common words.

In the plot below, the highest tf-idf words for each party are visualized, showcasing the terms that differentiate their plans from the others.

The TF-IDF analysis highlights the unique focus of each political party. Acción Popular emphasizes pension systems (e.g., SNP, SPP, CIC) and references historical figures like Belaúnde. Fuerza Popular focuses on public health and pandemic-related terms like infectados (infected) and defendemos (we defend), along with political events such as Noviembre (November protests). Juntos por el Perú centers on social rights and minorities, with words like LGBTI, Afroperuanos, and solidarias (solidarity).

Partido Morado highlights innovation and societal issues, referencing talento (talent), deporte (sports), and maltrato (abuse). Podemos prioritizes economic relief with terms like IGV (tax), deudas (debts), and planillas (payrolls). Victoria Nacional focuses on infrastructure and judicial reforms, with terms like Ernc (renewable energy), crudo (crude oil), and jurados (juries). These unique keywords underline each party’s strategic priorities.

There are some abbreviations:

Party	Acronym	Meaning
Acción Popular	SNP	National Pension System (Sistema Nacional de Pensiones)
	SPP	Private Pension System (Sistema Privado de Pensiones)
	OECD	Organization for Economic Cooperation and Development (Organización para la Cooperación y el Desarrollo Económicos)
	FCRP	Consolidated Reserve Fund of Pensions (Fondo Consolidado de Reservas Previsionales)
	CIC	Individual Capitalization Account (Cuenta Individual de Capitalización)
JP	APS	Primary Health Care (Atención Primaria de Salud)
	LGBTI	Lesbian, Gay, Bisexual, Transgender, and Intersex Community (Comunidad: Lesbianas, Gays, Bisexuales, Transgénero e Intersexuales)
	SSR	Sexual and Reproductive Health (Salud Sexual y Reproductiva)
	RDNFO	National Fiber Optic Backbone (Red Dorsal Nacional de Fibra Óptica)
	SPC	Peruvian Care System (Sistema Peruano de Cuidados - Proposal)
Podemos	FAE	Business Support Fund (Fondo de Apoyo Empresarial - FAE Agro)
	IGV	General Sales Tax (Impuesto General a las Ventas)
	LIC	Criminal Investigation Institute (Instituto de Investigación Criminal - Proposal)
	PYME	Small and Medium Enterprises (Pequeña y Mediana Empresa)
	EIRLS	Limited Liability Individual Enterprises (Empresas Individuales de Responsabilidad Limitada)
	Tele	Prefix (e.g., Tele-work or remote work)
Fuerza Popular	—	No specific acronyms mentioned
Partido Morado	—	No specific acronyms mentioned
Victoria Nacional	ERNC	Non-conventional Renewable Energies (Energías Renovables No Convencionales)
	CTI	Science, Technology, and Innovation (Ciencia, Tecnología e Innovación)
	Lesa	Crimes Against Humanity (Crímenes de Lesa Humanidad)
	FONDEM	Trust Fund for Financing Sustainable Balanced Development (Fondos de Fideicomiso para Financiar el Desarrollo Equilibrado Sostenible)
	—	Multisectoral Zones of Natural Resource Concessions (Multisectorial de Zonas de Influencia de Concesiones de Recursos Naturales - Proposal)

We can apply the same TF-IDF analysis using combinations of words (bigrams):

The bigram TF-IDF analysis highlights unique combinations of terms that define each party’s priorities. Acción Popular focuses on governance and opportunities with phrases like Fernando Belaúnde and Mejores Oportunidades. Fuerza Popular emphasizes educational and social programs, using terms like Programas Sociales and Brecha Digital (digital gap). Juntos por el Perú underscores social justice with bigrams like Derechos Colectivos (collective rights) and Diversidad Sexual (sexual diversity).

Partido Morado prioritizes innovation and sustainability, evident in phrases like Servicios Digitales (digital services) and Manifestaciones Culturales (cultural manifestations). Podemos highlights economic concerns with terms like Baja Productividad (low productivity) and Futuras Pandemias (future pandemics). Victoria Nacional focuses on infrastructure and legal reforms, with bigrams such as Gas Natural (natural gas) and Centros Poblados (populated centers). These combinations offer deeper insights into each party’s distinct focus areas.

So far, we have interpreted words as carriers of meaning, identifying those that are most relevant to each government plan. This has allowed us to highlight unique terms and ideas across the parties. Now, we will shift our focus to interpret words as carriers of sentiments and examine which government plans make greater use of positive or negative sentiments.

Sentiment Comparison

First, let’s identify the most frequently mentioned words associated with positive (Positivo) and negative (Negativo) sentiments in the government plans:

The plot reveals a mix of positive and negative sentiments across the parties’ government plans. Positive sentiments dominate with words like reforma (reform), derechos (rights), sostenible (sustainable), and capacidades (capacities), highlighting a focus on social change, sustainability, and empowerment. Negative sentiments, while less frequent overall, are concentrated in terms like pobreza (poverty), crisis (crisis), and deudas (debts), reflecting concerns about socio-economic challenges.

Overall, parties like Juntos por el Perú and Partido Morado emphasize positive framing, focusing on opportunities and rights, while Acción Popular and Podemos lean more on negative sentiments, addressing pressing issues like poverty and financial instability. These contrasts underscore the varied approaches to addressing Peru’s challenges.

This visualization highlights the words contributing most significantly to positive and negative sentiments in each government plan. The distribution of these words offers insight into the tone and emotional framing each party uses in their proposals.

By summing the number of positive words and subtracting the number of negative words in each government plan, we can calculate a score that provides an indication of the overall sentiment (positive/negative) expressed in each document.

From the results, it is evident that Partido Morado and Juntos por el Perú use more words that convey positive sentiments (e.g., derechos [rights], reforma [reform], sostenible [sustainable], integral [comprehensive], etc.). On the other hand, Podemos and Acción Popular use more words that express negative sentiments (e.g., pobreza [poverty], baja [decline], crisis [crisis], deuda [debt], etc.). This provides insight into the emotional tone each party aims to evoke through their government plan.

Correlation Between Parties

Finally, we can examine the correlation between the parties by analyzing the percentage of usage for each word and comparing how similar their usage patterns are across the government plans. If you cant see the page, click here.

This analysis reveals the likelihood that two government plans use similar words at comparable frequencies. By identifying these correlations, we can gauge how aligned or distinct the parties are in their language and focus areas. It reveals the linguistic similarities between the political parties’ government plans.

Victoria Nacional and Podemos share the highest correlation (0.8), suggesting significant overlap in their vocabulary and thematic focus. Partido Morado also shows high correlations with JP (0.8) and Podemos (0.7), indicating shared themes or priorities. Acción Popular exhibits the lowest correlations (0.6) across all comparisons, reflecting its distinct linguistic and thematic approach. These relationships highlight the varying degrees of alignment or divergence in the parties’ proposals.

Conclusions

Finally, we can draw comparisons between the different political parties based on their government plans. It is important to clarify that all political parties refer to general, demagogic topics such as Education, Health, Employment, Economic Recovery, etc. Therefore, the conclusions below aim to highlight the subtle characteristics that make each plan stand out.

Acción Popular

Unlike other government plans, Acción Popular has a notable emphasis on terms related to pensions (e.g., pension, SNP, SPP, CIC, affiliates, etc.), which is reflected in the frequency of these words.
Sentiment analysis indicates that the plan leans towards negative sentiments, focusing on words like poverty, debts, taxation, and abusive.

Fuerza Popular

Displays a strong presence of terms related to social assistance and support for small businesses (e.g., social, plan, programs, resources, etc.). Additionally, it highlights the current pandemic and the distribution of vaccines/tests.
Sentiment analysis shows a balanced use of positive and negative words. While it includes terms related to poverty (e.g., struggle, crisis, problems), it also incorporates a similar number of words related to support (e.g., rights, rescue, social programs, we defend).
Makes an uncommon reference (compared to other plans) to the November 2020 protests against Manuel Merino, emphasizing that the protests “negatively affected” the police’s public image.

Juntos por el Perú (JP)

Juntos por el Perú differs from other government plans by focusing heavily on social rights and assistance for minorities (e.g., social, cultural, communities, services, LGBTI, Afro-Peruvians, etc.).
Sentiment analysis indicates that the plan emphasizes positive sentiments, using terms such as rights, capabilities, and strengthen.
Includes topics related to education and sexual diversity, which are not widely addressed in other plans.

Partido Morado

This government plan closely resembles others by extensively addressing common topics in government plans (e.g., health, education, services, etc.) without adding or omitting much compared to other plans.
However, it also addresses unique topics, though in a less detailed manner, such as support for talent and sports, republicanism, leveraging the fourth industrial revolution, animal abuse prevention, support for families and adolescents, among others.
Sentiment analysis shows a strong emphasis on positive sentiments, making little reference (as a proportion of total words) to social problems or negative words such as poverty, risk, or diseases.

Podemos

The plan prominently references tax relief and business support (e.g., IGV, businesses, debts, payrolls, SMEs [PYME], limited liability enterprises [EIRLS], mortgages, collections, etc.). Additionally, it highlights remote services (e.g., telemedicine, tele-education).
Sentiment analysis indicates a strong emphasis on negative sentiments, frequently using words related to shortcomings and problems (e.g., decline, crisis, debts, etc.).

Victoria Nacional

This government plan places a strong focus on projects (e.g., gas, crude oil, supply chains, infrastructure, etc.) and judicial topics (e.g., law, juries, crimes against humanity, penal court, international criminal court, etc.). Notably, Forzyth has recently mentioned plans to reform the Peruvian judicial system.
Sentiment analysis shows a balanced use of positive and negative sentiments. While it includes words related to problems (e.g., lack, poverty), it also incorporates an equal number of terms related to projects (e.g., reform, sustainable, competencies, etc.).

Final Conclusion

This text analysis of the 2021 Peruvian presidential election government plans offers a detailed comparison of the linguistic and thematic elements used by the six major political parties. The primary goal was to uncover distinct patterns, highlight the emotional tone, and identify the focus areas of each party while remaining impartial.

Tools Used

The analysis was conducted using the following tools and libraries:

Text Mining Tools: tidytext, widyr, igraph, ggraph, and wordcloud2 for processing, tokenization, and visualization of textual data.
Visualization Libraries: ggplot2 and GGally for plotting sentiment trends, correlations, and word frequencies.
Sentiment Analysis: Using a custom sentiment dictionary, the analysis classified words into positive or negative categories and calculated sentiment scores for each party.
TF-IDF Analysis: To highlight words unique to each government plan, we implemented the Term Frequency-Inverse Document Frequency (TF-IDF) methodology, filtering out common terms across all plans.
Data Wrangling: tidyverse libraries for organizing, cleaning, and summarizing data.

Key Insights

Common Themes: Across all plans, terms related to education, health, economic recovery, and governance were prevalent, emphasizing their broad appeal and political salience.
Sentiment Distribution: While most parties addressed similar issues, their emotional framing varied. Juntos por el Perú and Partido Morado stood out for their positive sentiment, emphasizing rights and progressive reforms, whereas Podemos and Acción Popular leaned towards negative sentiments, focusing on crisis and problems. This divergence highlights strategic attempts to evoke specific emotional responses from their audience.
Distinct Focus Areas:
- Acción Popular emphasized pension reform and taxation.
- Fuerza Popular focused on social assistance and small business support, including pandemic-related issues.
- Juntos por el Perú highlighted minority rights and social equality, with a unique focus on education and sexual diversity.
- Partido Morado addressed common governance topics but also touched on innovation, republican values, and family support.
- Podemos prioritized tax relief and remote services like telemedicine and tele-education.
- Victoria Nacional focused on infrastructure projects and judicial reform.
TF-IDF Highlights: This analysis employed advanced text-mining techniques such as TF-IDF to identify distinct patterns in government plans, allowing for a more nuanced understanding beyond mere descriptive analysis. These techniques were essential for filtering out generic terms and spotlighting unique priorities, enabling the extraction of actionable insights. Additionally, sentiment analysis provided a structured way to interpret the emotional undertones in each document, offering a comprehensive picture of the strategies employed by political parties. Together, these approaches demonstrate the power of data-driven methods in uncovering meaningful patterns in complex textual data.
Correlations: Partido Morado showed the highest overlap in vocabulary with other parties, indicating thematic alignment, while others demonstrated distinct linguistic patterns.

Implications

This analysis provides voters and researchers with a structured overview of the linguistic and thematic strategies employed in political discourse. It highlights the importance of sentiment and word frequency analysis in understanding how political parties frame their priorities and connect with their audience. Additionally, by identifying unique terms through TF-IDF, we can better distinguish the distinctive focuses of each party, beyond the generic promises common in election campaigns.

These results demonstrate the potential to extract key terms from diverse texts, providing valuable input for Machine Learning or AI systems. By leveraging this structured data, such systems can organize and interpret multiple documents more accurately, minimizing the risk of hallucinations and enhancing overall reliability.

By combining qualitative interpretation with robust text analysis tools, this study underscores the power of data-driven approaches in political and social research.

Final Note: This analysis reveals a common bias in political research: an overemphasis on urban areas and mainstream perspectives. The victory of a smaller party campaigning in rural Peru demonstrates the importance of considering diverse regions and populations in political studies. This lesson emphasizes the need for inclusive research methods that capture the full complexity of political dynamics.

Election Results Update

The 2021 Peruvian general elections resulted in a significant political shift. Pedro Castillo, representing Perú Libre, narrowly won the presidency with 50.13% of the votes in the second round, defeating Keiko Fujimori of Fuerza Popular, who secured 49.87%. This tight margin highlights the country’s deep political polarization.

In the congressional elections, the distribution of seats presents a fragmented political landscape. Perú Libre emerged as the largest single party with 37 seats, followed by Fuerza Popular with 24 seats. Other notable parties include Acción Popular (16 seats), Alianza para el Progreso (15 seats), and Renovación Popular (13 seats). Smaller parties such as Juntos por el Perú, Somos Perú, and Podemos Perú captured fewer seats, reflecting diverse political representation.

This demonstrates the electorate’s divided support across a wide range of ideologies and priorities, showcasing the complexity of Peru’s political landscape. The fragmented composition of Peru’s Congress, with no party holding a clear majority, has significantly contributed to the nation’s political instability in recent years. This division has made coalition-building and legislative consensus challenging, reflecting the electorate’s diverse ideological preferences. The lack of a dominant party in Congress set the stage for subsequent turmoil, including President Pedro Castillo’s failed attempt to dissolve Congress in December 2022, which led to his impeachment and arrest. This event further exacerbated the political crisis, resulting in a government perceived as lacking broad representation and struggling to achieve stability.

Among the parties analyzed in the blog, Fuerza Popular, Acción Popular, Podemos Perú, Juntos por el Perú, and Partido Morado were central to the textual and sentiment analysis. Their performance in the election was varied, reflecting the disconnect between campaign messaging and voter behavior.

Fuerza Popular secured 24 seats in Congress, demonstrating its ability to maintain a significant political presence. Its campaign, which balanced both negative and positive sentiments, might have appealed to voters seeking stability amid economic and social challenges.

Fuerza Popular’s discourse reflects a balance between addressing immediate societal challenges and proposing support-based solutions. Key words like “pobreza” (poverty), “crecimiento” (growth), “rescate” (rescue), and “defendemos” (we defend) reveal a focus on societal problems and the promise of interventions. Additionally, terms such as “segmentos poblacionales” (population segments) and “programas sociales” (social programs) suggest a tailored approach to addressing diverse needs, potentially resonating with voters concerned about inclusivity and systemic support. The emphasis on “infectados” (infected) and “pruebas” (tests) highlights an active response to pandemic-related challenges, a pressing issue during the election. These targeted messages, combined with a balance of positive sentiments such as “derecho” (rights) and negative sentiments addressing critical issues, may have attracted voters seeking pragmatic and immediate solutions to national crises. This blend of problem acknowledgment and actionable solutions likely positioned Fuerza Popular as a relatable choice for many Peruvians.

Acción Popular, with 16 seats, highlighted pensions and social reforms, consistent with its focus on words like reforma (reform) and sostenible (sustainable). However, its emphasis on negative sentiments may have limited broader appeal.
Podemos Perú captured only five seats, possibly due to its strong focus on tax relief and economic challenges but a heavier reliance on negative sentiments like deudas (debts) and crisis (crisis).
Juntos por el Perú obtained five seats despite its strong emphasis on social rights and positive sentiments, which were reflected in its use of terms like derechos (rights) and capacidad (capacity). This outcome suggests that its appeal may have been concentrated in specific demographic or ideological niches.
Partido Morado, with three seats, showed alignment with themes of innovation and republicanism but struggled to resonate broadly, as indicated by its minimal representation in Congress.

The Partido Morado focused on a highly positive discourse, emphasizing concepts like “reforma” (reform), “sostenible” (sustainable), “servicio” (service), and “derechos” (rights), which aligned with progressive and optimistic sentiments. Its unique themes included “servicios digitales” (digital services), “manifestaciones culturales” (cultural manifestations), and “republicanismo” (republicanism), as well as topics like animal welfare and environmental concerns. These ideas reflected a forward-thinking and modern platform.

However, the party’s discourse may have struggled to resonate with the broader Peruvian electorate for a few reasons. First, its topics, while innovative, might have seemed distant or abstract to voters concerned with immediate and tangible issues like poverty and crisis, which were heavily addressed by other parties like Fuerza Popular. Second, the highly technical or niche focus on topics like “datos abiertos” (open data) and “biomasa” (biomass) might not have connected with the needs or priorities of rural or economically vulnerable populations. Lastly, the lack of significant negative sentiments in their plan may have made it seem disconnected from the reality of many voters who were experiencing hardship and expected a party to explicitly address such struggles.

Ultimately, Partido Morado’s progressive and technocratic approach, while valuable, might have been perceived as elitist or overly idealistic, limiting its appeal among a diverse electorate facing immediate socio-economic challenges.

These results underscore the complexities of Peruvian voter preferences, which favored Perú Libre—a party not initially included in the blog’s analysis. This demonstrates the challenges of relying solely on media and surveys centered in Lima, highlighting the importance of addressing the full spectrum of voter priorities, particularly in rural areas.

Hello R Markdown

Thu, 01 Dec 2022 21:13:14 -0500

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

You can embed an R code chunk like this:

summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
fit <- lm(dist ~ speed, data = cars)
fit
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Coefficients:
## (Intercept) speed
## -17.579 3.932

Including Plots

You can also embed plots. See Figure 1 for example:

par(mar = c(0, 1, 0, 1))
pie(
c(280, 60, 20),
c('Sky', 'Sunny side of pyramid', 'Shady side of pyramid'),
col = c('#0292D8', '#F7EA39', '#C4B632'),
init.angle = -50, border = NA
)

Figure 1: A fancy pie chart.

Webscraping 201

Sun, 19 Sep 2021 00:00:00 +0000

More and more central government institutions are centralizing their main websites on the www.gob.pe platform. While the government also has a dedicated data website (Open Data Platform), unstructured information (like reports or PDFs) often contains more details.

For this example, we will retrieve the links to daily temperature reports from the main ports of the country, which are published daily in the Daily Oceanographic Bulletin by IMARPE and made available on the government website. Since individual links do not follow a simple query structure, it is not possible to predict the daily URL. Therefore, we will use the government’s search tool to obtain the daily links.

Gob.pe Search Tool

To scrape information from the government website, we will use its search tool. This tool also operates with a query mechanism, simplifying the search process.

To review the parameters used, visit www.gob.pe and search for the desired information. At first glance, the query starts with the term /busquedas?.

For example, if we need to view IMARPE’s latest publications within a defined interval, we use the query:

We start with /busquedas? to indicate we are using the search tool, followed by the type of content, which in this case is publications: contenido[]=publicaciones, and the institution: institucion[]=imarpe.

Next, we specify the search parameters. To define a time interval, we use the parameters desde (from) and hasta (to): desde=04-07-2018&hasta=16-09-2021. Since IMARPE also publishes weekly reports, we include a term to filter results with the word “diario” (daily): term=Diario.

Thus, the data request URL would be as follows:

https://www.gob.pe/busquedas?contenido[]=publicaciones&institucion[]=imarpe&desde**=04-09-2021**&hasta**=16-09-2021**&term=**Diario**

Entering this address into a web browser directs us to IMARPE’s latest daily bulletin publications. This approach can be adapted to various types of information across institutions by experimenting with the search tool and reviewing the query.

This allows us to access the daily bulletin publications.

Webscraping with Rvest

With the web address defined, we proceed to retrieve the links to each publication.

path = "https://www.gob.pe/busquedas?contenido[]=publicaciones&institucion[]=imarpe&reason=sheet&sheet=1"

For web scraping, we will use the Rvest library, which provides tools to read the HTML content of a webpage and extract information, such as URLs. The tool to download HTML content from a web address is read_html.

library(tidyverse)
library(lubridate)
library(stringr)
library(rvest)
imarpe = read_html(path)
imarpe

## {html_document}
## <html lang="es-pe">
## [1] <head>\n<meta http-equiv="Content-Type" content="text/html; charset=UTF-8 ...
## [2] <body>\n<!-- Google Tag Manager (noscript) -->\n <noscript><iframe s ...

Using the read_html function on the Peruvian government results webpage, we obtain a set of HTML information stored in the variable imarpe.

The variable imarpe now contains a set of instructions for web browsers, organized in various tags. To extract the tags, we use the html_head function, and to retrieve the text within these tags, we use the html_text function. Next, we use text extraction functions (refer to Regex on Google) to extract the URLs of IMARPE publications using text patterns.

listas = imarpe %>%
html_node("head") %>% # Extract URLs from the head section
html_text() %>%
str_extract_all(paste("href(.*?)",year(today()),sep="")) %>% # Pattern: Starts with href and ends with 2021
as.data.frame(col.names = "Text") %>%
as_tibble() %>%
filter(str_detect(Text, pattern = "diario")) %>% # Filter URLs containing "diario"
mutate(Text = paste("https://www.gob.pe",
str_remove_all(Text, "href=\\\\\\\""),
sep = "")) %>% # Convert to full web address
mutate(Date = dmy(str_sub(Text, -10, -1))) %>% # Extract the date
arrange(Date)
imarpe %>%
html_node("head") %>% # Extract URLs from the head section
html_text() %>%
str_extract_all("href(.*?)2024") %>% # Pattern: Starts with href and ends with 2021
as.data.frame(col.names = "Text") %>%
as_tibble() %>%
filter(str_detect(Text, pattern = "diario"))

## # A tibble: 2 × 1
## Text
## <chr>
## 1 "href=\\\"/institucion/imarpe/informes-publicaciones/6204259-boletin-diario-o…
## 2 "href=\\\"/institucion/imarpe/informes-publicaciones/6200693-boletin-diario-o…

To process the text, we must first review the HTML text and its characteristics. The text processing functions are as follows:

First, extract text (str_extract_all) that starts with href (the tag indicating a URL in HTML) and ends with 2021 (the web addresses of interest include the year in their query). Using Regex, the dot (.) matches all characters. The asterisk (*) matches one or more, while the question mark ensures that the Regex is not greedy (captures the minimum number of characters possible to end before 2021).
Convert the extracted vectors into a dataset.
Filter URLs containing the word “diario” since we are looking for IMARPE’s Daily Oceanographic Bulletins, which include this word in their URLs.
Remove the href HTML tag from the text using the str_remove_all function.
Extract the dates by taking the last 10 characters of the URL (the query includes the date at the end of the URL).

listas

## # A tibble: 2 × 2
## Text Date
## <chr> <date>
## 1 https://www.gob.pe/institucion/imarpe/informes-publicaciones/62006… 2024-11-19
## 2 https://www.gob.pe/institucion/imarpe/informes-publicaciones/62042… 2024-11-20

This provides a list of URLs, which can later be used to download each publication.

Rvest is not the only method for web scraping with R. Additionally, more advanced tools, such as headless browsers (simulated web browsers capable of clicking and entering data like a normal browser), can be used. Examples include Rselenium (more complex) or Webdriver (simpler).

Conclusion

Using web scraping techniques and text manipulation, it is possible to extract the URLs of IMARPE publications via the search tool on the government website (www.gob.pe).

Using similar web scraping techniques, each publication can be accessed through a download loop to retrieve the uploaded information, which will be demonstrated in a future post.

Office Skills with Word

Fri, 02 Jul 2021 00:00:00 +0000

Tips for Office Skills in Microsoft Word

Word is the text processing tool in Microsoft Office, frequently used to create documents. While it is very user-friendly, there are some tips to accelerate text editing. This can be particularly important when regularly working with a large number of documents.

Word Shortcuts

You can select a sentence within a paragraph by pressing “Ctrl + Click”. This allows you to quickly select a sentence, for example, when you want the first sentence to be bold (Ctrl + B).
You can move any selection by dragging it with the mouse.
Use the combination “Alt + Shift + Arrow Up/Down” to move a paragraph.
Justify a paragraph with the combination “Ctrl + J.”
Adjust line spacing using different combinations: “Ctrl + 1” (single spacing), “Ctrl + 2” (double spacing), “Ctrl + 5” (1.5 spacing), and “Ctrl + 0” (adds a blank line before the paragraph).
Delete text faster (word by word instead of letter by letter) by pressing “Ctrl + DELETE” or “Ctrl + BACKSPACE”.
The Clipboard stores all copied items in chronological order. Access it by clicking the arrow in the lower-right corner of the Clipboard section under Home.
Insert a page break with the combination “Ctrl + Enter.”

Styles

Manually changing formatting repeatedly can be inefficient. You can save formatting to reuse in the future. Styles are located under the Home tab in the Styles section.

Modify a style by right-clicking it and selecting “Modify.” You can adjust indentation, justification, font type, font size, and more. You can also select all instances (paragraphs) where the style is used.

The following video explains the process well:

You can also update a style to match an existing paragraph. Select the paragraph and, in the Styles section, right-click the style, choosing “Update … to Match Selection”.

Additionally, you can create a new style. Select the paragraph with the modified style, and a window will appear with options. Click on Styles and then Create Style. This allows new text to adopt the predefined style.

Styles can be saved for the current document or for all documents created using the same template.

Locking Styles

An organization may require a corporate format template and prevent modifications to the formatting. In this case, you can lock the style settings. Note that this will restrict changes to font type, size, spacing, indentation, etc., except for the predefined styles. Ensure all required styles are created before locking.

To lock formatting, go to Review, click on Restrict Editing under Protect, and a window will open on the left. Choose Formatting Restrictions and its settings.

In the settings, deselect all styles except the corporate ones previously created.

Finally, you will see three format options, including “Block Theme or Scheme Changes” and “Prevent Quick Style Set Changes.” The first prevents theme changes (set of styles), and the second blocks changes to quick styles. It is recommended to enable both.

Press “Yes, Apply Protection” to finalize. Optionally, Word allows you to create a password to prevent disabling the formatting restrictions.

Ruler

Clicking on the ruler (above the document) will create an “L” marker. This marker sets where the cursor moves after pressing Tab. Clicking on the symbol to the left of the ruler changes it to an inverted “T.” Clicking on the ruler again adds another marker, indicating where the cursor will move (centered). This is a simple way to center titles and content. Experiment with various symbols and ruler functions.

Titles for Figures/Tables

You can add titles to figures and tables. The title formatting can also be modified under the Styles tab. This allows for numbering figures and tables and implementing a list of figures and tables.

Sources:

Webscraping 101

Sun, 28 Mar 2021 00:00:00 +0000

The internet offers a wealth of free and easily accessible information. There are two types of available data: data that can be easily downloaded through a “download” button and data visible on web pages.

To obtain both types of data, it is not necessary to manually download them. It is possible to automate the process through code.

In this post, I will introduce how to download the first type of data (data requiring only pressing a “download” button) automatically via code. However, to achieve this, it is important to understand the mechanisms that enable downloading data online. Today, we will discuss two such mechanisms: the POST and GET methods.

HTTP: GET and POST

The HTTP protocol is the communication protocol used by browsers to access web pages. HTTP regulates how the server (where the web page is hosted) sends resources to the client (the web browser). These resources contain the “instructions” a web browser (e.g., Chrome) uses to display and interact with the web page.

A URL is a web address that defines how a resource is located on the internet (e.g., https://www.google.com). When entering this address into a browser, it sends an HTTP request to the server, asking for the web resource (using a GET or POST method). The server receives the request, searches for the file in question (an HTML page, Excel file, etc.), and sends a header back to the browser. The header is essentially a message indicating whether the search was successful (whether the file exists). If successful, the server also sends the requested file.

For data downloads, two primary HTTP methods are commonly used: GET and POST. The GET method is used to “retrieve” resources from the server via a request. The POST method allows for “retrieving” and “saving” resources on the server by including additional data in the request.

GET: This HTTP method is used to retrieve information from a server. It sends a request to obtain a response (e.g., an Excel file or data we want to download) to the client (web page). A key feature of the GET method is that it allows the data to be transmitted “visibly” in the browser through the URL using queries. A query is a message added to the URL specifying the required information. For example, suppose we want to request visitor data for November from a hypothetical website, www.luis.com. The variables might include variable = visits and month = November. The query would look like: www.luis.com/download?variable=visits&month=november. The bold section is called the query. This query sends the request to the server and receives a response in return (e.g., an Excel file). Each website has its own syntax for how to structure such requests.
POST: Another HTTP method, which, unlike GET, sends information to be processed and stored on the server. The POST method can also be used to download data. The uniqueness of POST for downloading data lies in the use of two URLs. This improves server security. To download a file with the POST method, a query is sent similar to GET. However, instead of immediately providing a response, the server opens a temporary link at another URL for the file download.

Example of Downloading with GET and POST

To better understand the GET and POST protocols, we will use both methods to download data available on the COES website. The site contains data that can be downloaded using both methods.

Example: POST Method

The COES Indicators Portal allows for downloading daily electricity production data up to the previous day.

To access this section, go to the Indicators Portal on the COES website. Clicking the export button sends a query containing the selected date range in the request, and in return, an Excel file with the data is provided. How do we determine whether the method is POST or GET? It’s easy, using web analysis tools like the Chrome Developer Tools.

Figure 1: Chrome Developer Tools

Clicking on it reveals several tools, the most relevant being Network. Selecting Network shows the interaction between our client (web page) and the COES server. First, click Clear to clean up the workspace and focus on future interactions. Then, click the Export button on the COES webpage (after selecting the date range) to download the Excel file. In the Network area, two interactions appear. The presence of two interactions suggests the use of a POST method (one URL for sending and another for downloading).

Figure 2: Network

Click the first interaction (exportargeneracion) for more details. Since it’s the first interaction, it should contain the request sent to the server.

Figure 3: Headers - POST

We can see that the Request Method is POST. Additionally, the request URL is "https://www.coes.org.pe/Portal/portalinformacion/exportargeneracion". Finally, the request sends three pieces of data: the start date (fechaInicial), the end date (fechaFinal), and an indicator (indicador).

Note the date format: day/month/year, where both the day and month always have two digits. This format is crucial when sending the request.

A header acts as metadata (additional descriptive information) included in the request (request) and response. There are 15 request headers in the image. While not mandatory, they can be helpful, as we will see later.

Finally, the Status Code = 200 indicates a successful interaction.

Thus, we can conclude that the method for obtaining this information is POST, including the requested date range in the query.

Since it’s a POST method, the request does not provide an immediate response. We must check the second interaction to find the URL for downloading the data.

Clicking on the second interaction (descargargeneracion) reveals the following information:

Figure 4: Headers - GET

Here, a “gateway” to the data was opened at the URL: “https://www.coes.org.pe/Portal/portalinformacion/descargargeneracion”. This interaction uses the GET method, meaning that entering the URL automatically retrieves the data from the server. Since the request data was already sent in the previous interaction, no additional query is needed.

Example: GET Method

In this example, we will download data from the Daily Operation Evaluation Report (IEOD). This report provides more granular daily electricity consumption data, such as demand by geographic zones, large companies, resources, and more. To access this data, first, visit the IEOD platform and verify whether the download function uses POST or GET.

Figure 5: IEOD

After selecting a date on the IEOD platform and choosing the desired Excel file ("Anexo1_Resumen_0709.xlsx"), right-click on the link to copy the URL. In this case, the copied URL is: “https://www.coes.org.pe/portal/browser/download?url=Post%20Operaci%C3%B3n%2FReportes%2FIEOD%2F2020%2F09%20Setiembre%2F07%2FAnexo1_Resumen_0709.xlsx”. From the bold portion, it is evident that this is a GET method (data added to the request is visible in the URL, unlike POST). But which parts of the request change, and what do characters like %2F, %20, or %C3%B3 mean?

Using the Network section of the Chrome Developer Tools, we can see that the query starts after download?url=.

Figure 6: Headers - GET

From this, we can infer that the GET portion carrying the data, or the query, is: “Post Operación/Reportes/IEOD/**2020

Electricity for Economists I

Sun, 14 Mar 2021 00:00:00 +0000

Electricity demand is an excellent predictor of GDP. This information has a lag of only one day, which can help us quickly predict whether current GDP is growing or declining. However, it is essential to understand some basic concepts before using it.

Technical Aspects

Electric energy is defined as the movement of electrons through a conductor over a given period. The physical force or pressure that induces this movement is called voltage, measured in volts (V). The rate at which the electrons flow is called current intensity, measured in amperes (I).

Electric current (I) is produced by the movement of electrons caused by a voltage difference (V).

The electric power, measured in watts (W), quantifies the amount of energy consumed, produced, or transferred per unit of time. Electric energy, on the other hand, represents the total amount of energy consumed, produced, or transferred during a specific period and is measured in watt-hours (Wh).

For example, if the power of a lightbulb is 100 W and it remains on for two hours, the electric energy consumed would be 200 Wh.

Figure 1: Electricity

Figure 2: Energy vs Power

It’s common to find these units with the following prefixes. We will frequently use the giga unit when conducting analyses:

Prefix	Symbol	`$10^n$`	Equivalent
kilo	K	`$10^3$`	1,000
mega	M	`$10^6$`	1,000,000
giga	G	`$10^9$`	1,000,000,000

Demand

Electricity’s utility doesn’t derive from its direct consumption but from its ability to power equipment, making it a derived demand from other economic agents' needs (industries, households, and governments).

In Peru, there are two types of users in the national electricity system:

Regulated Users: These are agents with a demand of less than 200 KW, typically residential users and small businesses. Their prices are regulated by the government entity OSINERGMIN.
Free Users: These are agents with a demand of over 200 KW, usually large factories and mines that can negotiate electricity prices. Users with a demand between 200 KW and 2500 KW can choose to be either regulated or free.

Generation, Transmission, and Distribution

Electricity cannot be stored on a large scale at viable costs, necessitating its simultaneous production and consumption. This requires an infrastructure capable of generation, transportation, and distribution in real time.

Figure 3: Generation and Distribution

Generation

This first step transforms other energy forms (coal, gas, water flow) into electric energy. Generation relies on diverse technologies such as hydroelectric, thermal, solar, and wind plants, depending on the market size and energy resource availability.

In Peru, hydroelectric energy has historically dominated, but the development of the Camisea natural gas field has expanded gas-based electricity generation.

Figure 4: Generacion

Transmission

This step transports electricity from generation centers to consumption areas via high-capacity transmission lines. Due to economies of scale, this segment exhibits natural monopoly characteristics.

Figure 5: Cables de Transmision

In Peru, most transmission lines are 138 kV, 220 kV, and 500 kV. Higher capacity lines allow for greater energy transport, which is crucial for supporting industrial growth.

Distribution

Distribution involves delivering electricity to end-users through medium- and low-voltage networks. In Peru, this segment is often managed by state companies, except in large cities like Lima, where private companies like Luz del Sur operate.

Key Institutions

Ministry of Energy and Mines: Responsible for legislating energy-related matters.
OSINERGMIN: Supervises compliance and sets regulated user prices.
COES (Committee for Economic Operation of the System): Coordinates the electricity system, manages the short-term market, and plans long-term energy needs.

Figure 6: Modelo del Mercado Peruano

The COES exists to maintain the proper balance between energy demand and supply, managing the agents involved (Generators-Distributors) through the electricity market.

The decisions of the COES are mandatory for its members (generation, transmission, and distribution companies). Electricity generation is dispatched in order of the marginal cost of production. As a result, the first energy to be dispatched comes from the most efficient technologies (solar, wind, etc.).

Figure 7: COES Despacho

Electricity for Economists II

Sun, 14 Mar 2021 00:00:00 +0000

Sources of Information in Peru

COES

The COES manages the technical aspect of electricity flow (supply and demand) in the National Interconnected Electric System or SEIN. If you want to review the technical aspects of electricity, you can find them here.

As mentioned earlier, COES manages information on electricity generation, transmission, and demand in real time, ensuring the electrical system remains balanced (Supply = Demand).

COES also publishes this information on its website https://www.coes.org.pe/portal/.

There is a variety of information available on the COES website. Below, we will describe the most relevant aspects for our work:

Indicators of the Peruvian Energy Sector

This is the most accessible information on the COES website. There is a link on the homepage:

Figure 1: Indicador de generación total

This is the most accessible information on the COES website. There is a link on the homepage:

Figure 2: Total Generation Indicator

Here, we can access various information provided by agents in the electricity sector, such as:

Electricity generation by all Peruvian generating companies within the SEIN.
An estimate of electricity demand.
Announcements of failures and repairs within SEIN.
Information on the volume and capacity of hydroelectric company reservoirs.

The most relevant information for day-to-day economic analysis is electricity generation. Users can select a date range and view data on the website or download it as an Excel file. The website provides executed measurements every 30 minutes and projections for the rest of the current day. For dates prior to today, only executed data is available.

Figure 3: View Generation Data

Information is categorized by generating company and resource type (hydroelectric, gas, etc.).

Operations

Within the options available in the website’s top menu, we find various operation-related features of COES.

This section provides access to COES operation data, including the most important options::

Database

Under the “Database” section, users can access automated data before it is published in official reports. The most relevant option for us is the Generation Meters:

Generation Meters: Contains information on monthly electricity generation in the country, based on the selected date range. The most important data is Total Energy (MWh). The data for the previous month is uploaded on the third day of the month.

This section provides official data on the previous month’s electricity generation directly and before publication elsewhere.

Operation Program

The “Operation Program” section includes projections for future electricity demand and corresponding generation planning. The projections are more reliable when closer to the current date. For instance, daily projections are highly reliable, whereas weekly and medium-term projections are less so.

This information is not very relevant to our day-to-day work and does not need to be mastered; however, it can be useful to know for certain special requests.

Electricity demand can be mathematically divided into two parts: an almost invariant trend and a variable component that fluctuates according to the economic cycle. It is generally considered that the consumption of free users (industries, mines, etc.) constitutes the variable component, while the consumption of regulated users (households) represents the invariant trend. For medium-term projections, both free-user and regulated-user demand are calculated. The Vegetative Demand is the term commonly used to describe the trend component of electricity demand.

To project demand, COES calculates vegetative demand using an ARIMA model, while the variable component is obtained through a survey of all free users. This combined approach allows for medium-term projections.

To view this projection, click on Medium-Term Operation Program, select the date, and in the 2_Reports section, download the file PMPO_“Date”.zip. Within the zip file, extract another zip file Demand_“Date”.zip and open the Excel file **Medium-Term Demand Projection_ARIMA_“Date”.xlsx**. Inside, you will find the projection. However, it is important to note that the free-user data is only obtained through an annual survey, making the projection less reliable. Nonetheless, it is useful for planning investments in the electricity sector.

International Exchanges

The international exchanges section contains information on the import and export of electricity. The most relevant option within this section is Energy Meter Logs, which provides data on electricity imports and exports measured in MWh.

Electricity is transmitted in loads, which can be either positive or negative depending on the direction, i.e., whether the demand originates from one side or the other. This allows the load direction to be measured at substations located on the border.

Peru’s most significant international electricity exchange is with the neighboring country of Ecuador. This exchange is closely linked to the energy conditions in the northern region of Peru (demand or surplus energy), which is usually related to its economic activity.

Single-Line Diagram

Finally, the most relevant feature within the operations section is the single-line diagram. The single-line diagram is a map (in PDF format) that contains the name and location of each electricity measurement point within the SEIN.

The map allows us to observe electricity transmission lines, which are subdivided into smaller lines feeding smaller sectors. This way, if we have more granular information on regionally focused demand (Post-Operation - Agent Demand - Distributors), we can determine with certainty whether one measurement point includes the data of another point (i.e., if it is the “parent” of another line), avoiding overestimating demand when aggregating all information.

Post-Operation

Within the options available in the top menu of the website, we find various features, including the post-operation options of COES. This section contains information processed after passing through the Operation section. It allows access to more granular and processed data on national electricity demand.

Figure 4: Post Operation

IEOD

The most important information we work with comes from the Daily Operation Evaluation Report (IEOD). This report contains all the detailed data from COES operations conducted on a specific day.

By accessing this section, we can select the desired date and download the most relevant Excel file: Anexo1_Resumen_“Date”.xlsx.

Part of the information divides the country into three regions: North, Center, and South, as shown in the following map:

Figure 5: Diagrama unifilar

In this context, the central region is defined as the area excluding the northern and southern regions. Iquitos is not part of SEIN and is considered an independent system.

Within the IEOD Excel file, we find various useful sheets, the most relevant of which are:

Demanda_UL: This sheet shows the power demand of Large Free Users (unlike free users [>2,500 KW], large free users are those demanding even more electricity [>10,000 KW]). The companies are categorized by their location in the North, Center, or South. Since the information is in power and provided every 30 minutes, to obtain the daily demand in GWh, the data must be summed, multiplied by 1,000, and divided by 2 (to convert the unit to 1 hour).

This data allows us to infer the demand of mining and manufacturing sectors. Companies are not legally obligated to provide this data (though there are monetary penalties), so the information may have gaps, represented by a demand of 0. Researchers must use their discretion to impute the missing data.

Demanda_Areas: This sheet shows the total power demand in the country, separated by Areas and Sub-Areas. The most reliable data is by Areas (North, Center, and South). To calculate the daily data, the values recorded every 30 minutes are summed and then divided by 2. The total corresponds to the energy observed by COES-affiliated companies, so it is slightly less than the total demand measured in SEIN (Indicators of the Peruvian Energy Sector).

These two sheets are the most relevant. However, the Excel file also includes information on generation by resource (hydroelectric, natural gas, etc.), which can be useful depending on the researcher’s needs.

Reports - Evaluation

This section displays the reports prepared by COES for the public. The most relevant report is the Monthly Report, which provides the official data on generation and the percentage variation for a given month. The monthly data does not match the sum of the daily data collected due to measurement errors.

The report is published on the 15th of each month. However, the official monthly generation data can be accessed earlier in the Operations -> Database section.

Agent Demand

This section shows the demand of two COES agents: Free Users and Distribution Companies:

Free Users: Displays the demand of some free users affiliated with COES. Each month, some free users voluntarily provide their daily demand data to COES, while others opt out, making this information variable. Compared to the IEOD, this section includes data from more companies and is subdivided by each company’s measurement points (e.g., each factory).
Distribution Companies: Shows the electrical demand in the region assigned to each distribution company. A measurement point may include data from another measurement point, which is why it’s important to review the Single-Line Diagram to ensure the information is consistent. The most relevant data in this section comes from higher-level measurement points (e.g., substations in the regional capital). To locate these, it is also useful to consult the single-line diagram. This information is somewhat complex to process and relatively new, so I conducted an analysis here.

Distribution companies are often state-owned and inefficient, so the data is variable and unreliable. This information is still being implemented and is not available for all districts. However, COES is working to ensure that, in the future, each distribution company will automatically report its demand data.

Members

Finally, another option in the top menu is “Members.” By clicking on “List of Members”, we can view the members of COES, whether they are Generators, Distributors, or Free Users.

Additionally, COES responds to any inquiries we may have. This is useful, for example, if we suspect a data error. To contact them, simply click on the blue button located at the top, labeled “Contact Us” when you hover over it.

OSINERGMIN

OSINERGMIN is the entity responsible for overseeing the members of COES. An interesting feature is the SEIN map, which contains the geographic location of SEIN components: https://www.osinergmin.gob.pe/newweb/uploads/Publico/MapaSEIN/.

Additionally, it is possible to access information on the demand of all free users in the country through the following link, under Publications -> Monthly Reports -> Free Electricity Market:

https://www.osinergmin.gob.pe/seccion/institucional/regulacion-tarifaria/publicaciones/regulacion-tarifaria

However, this information is delayed by 2 months, making it less useful for recent data needs.

MINISTRY OF ENERGY AND MINES

Finally, the Ministry of Energy and Mines also conducts an annual review of the electricity sector. Their reports are delayed by 1 or 2 years but contain a wealth of historical information, such as the total demand in the country by industrial activity.

The information can be downloaded from: http://www.minem.gob.pe/_estadisticaSector.php?idSector=6

Writing technical content in Academic

Fri, 12 Jul 2019 00:00:00 +0000

Academic is designed to give technical content creators a seamless experience. You can focus on the content and Academic handles the rest.

Highlight your code snippets, take notes on math classes, and draw diagrams from textual representation.

On this page, you’ll find some examples of the types of technical content that can be rendered with Academic.

Examples

Code

Academic supports a Markdown extension for highlighting code syntax. You can enable this feature by toggling the highlight option in your config/_default/params.toml file.

```python
import pandas as pd
data = pd.read_csv("data.csv")
data.head()
```

renders as

import pandas as pd
data = pd.read_csv("data.csv")
data.head()

Charts

Academic supports the popular Plotly chart format.

Save your Plotly JSON in your page folder, for example chart.json, and then add the {{< chart data="chart" >}} shortcode where you would like the chart to appear.

Demo:

You might also find the Plotly JSON Editor useful.

Math

Academic supports a Markdown extension for $\LaTeX$ math. You can enable this feature by toggling the math option in your config/_default/params.toml file.

To render inline or block math, wrap your LaTeX math with $...$ or $$...$$, respectively.

Example math block:

$$\gamma_{n} = \frac{
\left | \left (\mathbf x_{n} - \mathbf x_{n-1} \right )^T
\left [\nabla F (\mathbf x_{n}) - \nabla F (\mathbf x_{n-1}) \right ] \right |}
{\left \|\nabla F(\mathbf{x}_{n}) - \nabla F(\mathbf{x}_{n-1}) \right \|^2}$$

renders as

$$\gamma_{n} = \frac{ \left | \left (\mathbf x_{n} - \mathbf x_{n-1} \right )^T \left [\nabla F (\mathbf x_{n}) - \nabla F (\mathbf x_{n-1}) \right ] \right |}{\left |\nabla F(\mathbf{x}_{n}) - \nabla F(\mathbf{x}_{n-1}) \right |^2}$$

Example inline math $\nabla F(\mathbf{x}_{n})$ renders as $\nabla F(\mathbf{x}_{n})$.

Example multi-line math using the \\\\ math linebreak:

$$f(k;p_{0}^{*}) = \begin{cases}p_{0}^{*} & \text{if }k=1, \\\\
1-p_{0}^{*} & \text{if }k=0.\end{cases}$$

renders as

$$f(k;p_{0}^{*}) = \begin{cases}p_{0}^{*} & \text{if }k=1, \\
1-p_{0}^{*} & \text{if }k=0.\end{cases}$$

Diagrams

Academic supports a Markdown extension for diagrams. You can enable this feature by toggling the diagram option in your config/_default/params.toml file or by adding diagram: true to your page front matter.

An example flowchart:

```mermaid
graph TD
A[Hard] -->|Text| B(Round)
B --> C{Decision}
C -->|One| D[Result 1]
C -->|Two| E[Result 2]
```

renders as

graph TD
A[Hard] -->|Text| B(Round)
B --> C{Decision}
C -->|One| D[Result 1]
C -->|Two| E[Result 2]

An example sequence diagram:

```mermaid
sequenceDiagram
Alice->>John: Hello John, how are you?
loop Healthcheck
John->>John: Fight against hypochondria
end
Note right of John: Rational thoughts!
John-->>Alice: Great!
John->>Bob: How about you?
Bob-->>John: Jolly good!
```

renders as

sequenceDiagram
Alice->>John: Hello John, how are you?
loop Healthcheck
John->>John: Fight against hypochondria
end
Note right of John: Rational thoughts!
John-->>Alice: Great!
John->>Bob: How about you?
Bob-->>John: Jolly good!

An example Gantt diagram:

```mermaid
gantt
section Section
Completed :done, des1, 2014-01-06,2014-01-08
Active :active, des2, 2014-01-07, 3d
Parallel 1 : des3, after des1, 1d
Parallel 2 : des4, after des1, 1d
Parallel 3 : des5, after des3, 1d
Parallel 4 : des6, after des4, 1d
```

renders as

gantt
section Section
Completed :done, des1, 2014-01-06,2014-01-08
Active :active, des2, 2014-01-07, 3d
Parallel 1 : des3, after des1, 1d
Parallel 2 : des4, after des1, 1d
Parallel 3 : des5, after des3, 1d
Parallel 4 : des6, after des4, 1d

An example class diagram:

```mermaid
classDiagram
Class01 <|-- AveryLongClass : Cool
<<interface>> Class01
Class09 --> C2 : Where am i?
Class09 --* C3
Class09 --|> Class07
Class07 : equals()
Class07 : Object[] elementData
Class01 : size()
Class01 : int chimp
Class01 : int gorilla
class Class10 {
<<service>>
int id
size()
}
```

renders as

classDiagram
Class01 <|-- AveryLongClass : Cool
<<interface>> Class01
Class09 --> C2 : Where am i?
Class09 --* C3
Class09 --|> Class07
Class07 : equals()
Class07 : Object[] elementData
Class01 : size()
Class01 : int chimp
Class01 : int gorilla
class Class10 {
<<service>>
int id
size()
}

An example state diagram:

```mermaid
stateDiagram
[*] --> Still
Still --> [*]
Still --> Moving
Moving --> Still
Moving --> Crash
Crash --> [*]
```

renders as

stateDiagram
[*] --> Still
Still --> [*]
Still --> Moving
Moving --> Still
Moving --> Crash
Crash --> [*]

Todo lists

You can even write your todo lists in Academic too:

- [x] Write math example
- [x] Write diagram example
- [ ] Do something else

renders as

Write math example
Write diagram example
Do something else

Tables

Represent your data in tables:

| First Header | Second Header |
| ------------- | ------------- |
| Content Cell | Content Cell |
| Content Cell | Content Cell |

renders as

First Header	Second Header
Content Cell	Content Cell
Content Cell	Content Cell

Callouts

Academic supports a shortcode for callouts, also referred to as asides, hints, or alerts. By wrapping a paragraph in {{% callout note %}} ... {{% /callout %}}, it will render as an aside.

{{% callout note %}}
A Markdown aside is useful for displaying notices, hints, or definitions to your readers.
{{% /callout %}}

renders as

A Markdown aside is useful for displaying notices, hints, or definitions to your readers.

Spoilers

Add a spoiler to a page to reveal text, such as an answer to a question, after a button is clicked.

{{< spoiler text="Click to view the spoiler" >}}
You found me!
{{< /spoiler >}}

renders as

Click to view the spoiler

You found me!

Icons

Academic enables you to use a wide range of icons from Font Awesome and Academicons in addition to emojis.

Here are some examples using the icon shortcode to render icons:

{{< icon name="terminal" pack="fas" >}} Terminal
{{< icon name="python" pack="fab" >}} Python
{{< icon name="r-project" pack="fab" >}} R

renders as

Terminal
Python
R

An example preprint / working paper

Sun, 07 Apr 2019 00:00:00 +0000

Create your slides in Markdown - click the Slides button to check out the example.

Supplementary notes can be added here, including code, math, and images.

Display Jupyter Notebooks with Academic

Tue, 05 Feb 2019 00:00:00 +0000

from IPython.core.display import Image
Image('https://www.python.org/static/community_logos/python-logo-master-v3-TM-flattened.png')

print("Welcome to Academic!")

Welcome to Academic!

Install Python and JupyterLab

Install Anaconda which includes Python 3 and JupyterLab.

Alternatively, install JupyterLab with pip3 install jupyterlab.

Create or upload a Jupyter notebook

Run the following commands in your Terminal, substituting <MY-WEBSITE-FOLDER> and <SHORT-POST-TITLE> with the file path to your Academic website folder and a short title for your blog post (use hyphens instead of spaces), respectively:

mkdir -p <MY-WEBSITE-FOLDER>/content/post/<SHORT-POST-TITLE>/
cd <MY-WEBSITE-FOLDER>/content/post/<SHORT-POST-TITLE>/
jupyter lab index.ipynb

The jupyter command above will launch the JupyterLab editor, allowing us to add Academic metadata and write the content.

Edit your post metadata

The first cell of your Jupter notebook will contain your post metadata (front matter).

In Jupter, choose Markdown as the type of the first cell and wrap your Academic metadata in three dashes, indicating that it is YAML front matter:

---
title: My post's title
date: 2019-09-01
# Put any other Academic metadata here...
---

Edit the metadata of your post, using the documentation as a guide to the available options.

To set a featured image, place an image named featured into your post’s folder.

For other tips, such as using math, see the guide on writing content with Academic.

Convert notebook to Markdown

jupyter nbconvert index.ipynb --to markdown --NbConvertApp.output_files_dir=.

Example

This post was created with Jupyter. The orginal files can be found at https://github.com/gcushen/hugo-academic/tree/master/exampleSite/content/post/jupyter

Slides

Tue, 05 Feb 2019 00:00:00 +0000

Create slides in Markdown with Wowchemy

Wowchemy | Documentation

Features

Efficiently write slides in Markdown
3-in-1: Create, Present, and Publish your slides
Supports speaker notes
Mobile friendly slides

Controls

Next: Right Arrow or Space
Previous: Left Arrow
Start: Home
Finish: End
Overview: Esc
Speaker notes: S
Fullscreen: F
Zoom: Alt + Click
PDF Export: E

Code Highlighting

Inline code: variable

Code block:

porridge = "blueberry"
if porridge == "blueberry":
print("Eating...")

Math

In-line math: $x + y = z$

Block math:

$$ f\left( x \right) = ;\frac{{2\left( {x + 4} \right)\left( {x - 4} \right)}}{{\left( {x + 4} \right)\left( {x + 1} \right)}} $$

Fragments

Make content appear incrementally

{{% fragment %}} One {{% /fragment %}}
{{% fragment %}} **Two** {{% /fragment %}}
{{% fragment %}} Three {{% /fragment %}}

Press Space to play!

One Two Three

A fragment can accept two optional parameters:

class: use a custom style (requires definition in custom CSS)
weight: sets the order in which a fragment appears

Speaker Notes

Add speaker notes to your presentation

{{% speaker_note %}}
- Only the speaker can read these notes
- Press `S` key to view
{{% /speaker_note %}}

Press the S key to view the speaker notes!

Themes

black: Black background, white text, blue links (default)
white: White background, black text, blue links
league: Gray background, white text, blue links
beige: Beige background, dark text, brown links
sky: Blue background, thin dark text, blue links

night: Black background, thick white text, orange links
serif: Cappuccino background, gray text, brown links
simple: White background, black text, blue links
solarized: Cream-colored background, dark green text, blue links

Custom Slide

Customize the slide style and background

{{< slide background-image="/media/boards.jpg" >}}
{{< slide background-color="#0000FF" >}}
{{< slide class="my-style" >}}

Custom CSS Example

Let’s make headers navy colored.

Create assets/css/reveal_custom.css with:

.reveal section h1,
.reveal section h2,
.reveal section h3 {
color: navy;
}

Questions?

Ask

Documentation

An example journal article

Tue, 01 Sep 2015 00:00:00 +0000

Click the Cite button above to demo the feature to enable visitors to import publication metadata into their reference management software.

Create your slides in Markdown - click the Slides button to check out the example.

Supplementary notes can be added here, including code, math, and images.

An example conference paper

Mon, 01 Jul 2013 00:00:00 +0000

Click the Cite button above to demo the feature to enable visitors to import publication metadata into their reference management software.

Create your slides in Markdown - click the Slides button to check out the example.

Supplementary notes can be added here, including code, math, and images.

Mon, 01 Jan 0001 00:00:00 +0000

Electricity for for Economists III

Mon, 01 Jan 0001 00:00:00 +0000

Finally, we arrive at the third part of this series of blog posts introducing the electricity sector for economists. This time, we will explore how to automate the downloading and processing of electricity data to use it in our prediction models. As mentioned in previous posts, electricity information is useful due to its high correlation with GDP and its quick and easy accessibility.

Processing

The processing of COES data is done through web scraping: a technique used to extract data from web pages using software. Various techniques and software can be employed, but we will use the R software and HTML code to download the data.

If you want to review the information available on the COES website, you can do so here.

HTTP: GET and POST

To download information from the COES website (https://www.coes.org.pe/portal/), we will use HTML code. HTML is a markup language used for creating web pages. It contains the “instructions” that a web browser (such as Chrome) uses to display and interact with the web page.

For data downloads, two primary methods of data transfer in HTTP are commonly used: GET and POST.

GET: A method in the HTTP protocol used to retrieve information from a server. It sends a request to obtain a response (e.g., an Excel file or other data we want to download) for the client (web page). GET transmits the data in a “visible” way through the web browser’s URL. For example, suppose we want to request information from a hypothetical website, www.luis.com, and the required data for the request are the file name and month: name = luis and month = November. Using GET, the request would look like: www.luis.com/download?name=luis&month=november. Here, the bold part is referred to as the query. By submitting this request to the server (with our data), we would receive a response (e.g., an Excel file). Each website has its own syntax for submitting requests.
POST: Another method in the HTTP protocol that, unlike GET, sends information to the server for processing and updates the server. For example, when we submit (request) data through a form, the server processes the data and, after a redirect, returns (response) a page with information. To download a file using the POST method, a query is sent in a manner similar to GET. However, instead of immediately receiving a response, the server opens a temporary link at another URL for the file download.

Next, we will see a practical application of these terms to download relevant information from the COES website.

COES: Indicators of the Peruvian Energy Sector

To access this section, we go to the Indicators Portal on the COES website. By clicking on export, we send a query containing the selected date range in the request, receiving in return a response with an Excel file containing the data. How can we determine if the method is POST or GET? It’s easy, by using a web analysis tool such as the Chrome Developer Tools.

Figure 1: Figure 1

Clicking on the tools will display various options, with Network being the most relevant for us right now. By clicking on Network, we can observe the interactions between our client (web page) and the COES server. To clean up the workspace, first click Clear, which clears the window and helps focus on future interactions. Next, click on the Export button on the COES webpage (after selecting the date range) to download the Excel file. In the Network area, two interactions appear: having two interactions suggests that the method might be POST (one URL for sending and another for downloading).

Figure 2: Figure 2

Click on the first interaction (exportargeneracion) to gather more information. Since it’s the first interaction, it should contain the request sent to the server.

Figure 3: Figure 3

We can see that the Request Method is indeed POST. Additionally, the request URL is shown as “https://www.coes.org.pe/Portal/portalinformacion/exportargeneracion". Finally, we see that three data elements are sent in the request: the start date (fechaInicial), the end date (fechaFinal), and an indicator (indicador).

It’s important to note the date format: day/month/year, where both the day and the month always have two digits. This format is crucial when sending the request.

A header acts as metadata (additional descriptive information, though not intrinsic to the request) that accompanies the request. In the image, there are 15 request headers. While it’s not mandatory to use them, they can be helpful, as we will see later.

Finally, the Status Code = 200 indicates that the interaction was successful.

Thus, we can conclude that the method used to obtain this information was POST, with the requested date range included in the request.

However, since it’s a POST method, the request will not immediately provide a response. We need to check the second interaction to find out from which URL to download the data.

To do so, click on the second interaction (descargargeneracion) to access the following information:

Figure 4: Figure 4

Here, we can see that an “access gateway” to the data was opened at the URL: “https://www.coes.org.pe/Portal/portalinformacion/descargargeneracion". This interaction uses a GET method, meaning that by entering the URL, the server will automatically provide the data. Since the request data was already entered in the previous interaction, no additional data or query is needed.

Downloading in R

After analyzing the steps for downloading the files, we will now automate this process using R code. This will allow us to access the data directly through our code, enabling us to process it afterward.

First, R requires a special package to handle HTTPS methods, called httr. This package provides access to two essential functions, GET and POST, for downloading data. To install it, simply type install.packages(“httr”) once in the RStudio console.

After installing it, we invoke the library and use the POST method to download the file. Remember that the POST method requires two steps for file downloads: first, a POST request to send the data, and then a GET request to download the file from the temporary link that was opened. To download and read the data using the GET method, I will use the read.xlsx function from the openxlsx library, which allows downloading from a URL and reading an Excel file if contained within.

Since the link is temporary, the two functions must be executed sequentially. Below, I will explain these steps:

library(httr)
library(openxlsx)
fechaInicial = "07/09/2020"
fechaFinal = "07/09/2020"
POST(url = as.character("https://www.coes.org.pe/Portal/portalinformacion/exportargeneracion"), config(ssl_verifypeer = FALSE), add_headers("Referer"= "http://www.coes.org.pe/Portal/portalinformacion/generacion"),
body = list( "fechaInicial" = fechaInicial,
"fechaFinal"= fechaFinal,
"indicador"= 0)) ##Función POST envío los datos

## Response [https://www.coes.org.pe/Portal/portalinformacion/exportargeneracion]
## Date: 2024-11-23 19:28
## Status: 200
## Content-Type: application/json; charset=utf-8
## Size: 1 B

datos = read.xlsx("https://www.coes.org.pe/Portal/portalinformacion/descargargeneracion",sheet = 1,detectDates = TRUE) ##Funcion GET descargo los datos
datos[1:5,] #ejemplo de 5 primeras filas

## FECHA EMPRESA GRUPO/CENTRAL TIPO.GENERACIÓN
## 1 07/09/2020 AGRO INDUSTRIAL PARAMONGA PARAMONGA TERMOELÉCTRICA
## 2 07/09/2020 AGROAURORA S.A.C. MAPLE TERMOELÉCTRICA
## 3 07/09/2020 AGROINDUSTRIAS SAN JACINTO S.A.A. SAN_JACINTO TERMOELÉCTRICA
## 4 07/09/2020 AGUAS Y ENERGIA PERU PIAS1 HIDROELÉCTRICA
## 5 07/09/2020 ANDEAN POWER S.A.C. CARHUAC HIDROELÉCTRICA
## TOTAL
## 1 362.45
## 2 420.40
## 3 191.80
## 4 107.65
## 5 271.95

POST: The function uses the POST method, sending a request to the URL “https://www.coes.org.pe/Portal/portalinformacion/exportargeneracion", requesting the database we need between the dates fechaInicial (“07/09/2020”) and fechaFinal (“07/09/2020”). The query can be added in the body section as a list since there are multiple parameters. Finally, the indicator (“0”) seems to be a formality without any functionality. Note that we include the header “Referer”.

If the response in the console displays Status: 200, the request was successful. If this number does not appear, it means there was an error, and you should check the Chrome Developer Tools to ensure the parameters have not changed (parameters often change a couple of times per year).

The following are additional parameters, which are not mandatory but can help avoid errors:

A header acts as metadata (additional descriptive information, though not intrinsic to the request) included in the request. In this case, the header Referer refers to the page from which we redirected to the specified URL. If we were using a web browser, we would make the request from “http://www.coes.org.pe/Portal/portalinformacion/generacion” (Indicators Portal - Generation). While it is not necessary to include this, it informs the server that we accessed the data through a legitimate entry point and are not hackers attempting to overload their website or access unauthorized pages. Some websites have stricter security measures, while others are more lenient.

Finally, SSL Verifier refers to whether the website has anti-scam protection (e.g., preventing credit card numbers from being stolen). It is not mandatory to tell the function to ignore this, but we explicitly set it to FALSE just in case.

read.xlsx: The first parameter specifies the URL containing our Excel file. Since we performed the POST function earlier, the URL “https://www.coes.org.pe/Portal/portalinformacion/descargargeneracion" opened a temporary link where the requested Excel file is located. Sheet = 1 indicates which sheet of the Excel file to read. Finally, detectDates = TRUE ensures that Excel date formats are converted into standard date formats.

Processing in R

I will use the Tidyverse library, which provides several tools for data manipulation.

library(tidyverse)

To get an overview of the data:

glimpse(datos)

## Rows: 111
## Columns: 5
## $ FECHA <chr> "07/09/2020", "07/09/2020", "07/09/2020", "07/09/2020"…
## $ EMPRESA <chr> "AGRO INDUSTRIAL PARAMONGA", "AGROAURORA S.A.C.", "AGR…
## $ `GRUPO/CENTRAL` <chr> "PARAMONGA", "MAPLE", "SAN_JACINTO", "PIAS1", "CARHUAC…
## $ TIPO.GENERACIÓN <chr> "TERMOELÉCTRICA", "TERMOELÉCTRICA", "TERMOELÉCTRICA", …
## $ TOTAL <dbl> 362.45, 420.40, 191.80, 107.65, 271.95, 6.15, 242.40, …

The downloaded data contains five columns: Fecha (Date), Empresa (Company), Grupo (Group), Tipo de Generación (Type of Generation), and Total (MWh).

Typically, we only need to retain the Fecha and Total columns. However, the additional information can be useful for analysts (e.g., identifying generation from agricultural companies to infer sector activity or comparing gas, solar, and other types of generation).

Additionally, we might have downloaded data for a date range greater than one day. To aggregate total generation by day, we group the data by FECHA. Finally, we sum the values by day and convert MWh to GWh by dividing by 1,000.

datos %>% # Access the data
select(FECHA, TOTAL) %>% # Select Date and Total (MWh)
group_by(FECHA) %>% # Group by day
summarise(GWh = sum(TOTAL, na.rm = FALSE) / 1000) %>% # Sum and convert to GWh
ungroup() # Ungroup (optional)

## # A tibble: 1 × 2
## FECHA GWh
## <chr> <dbl>
## 1 07/09/2020 141.

This gives us daily data in GWh. The next step would be adding the data to a historical database. We won’t delve into that here but will instead move on to the next topic: IEOD.

COES: IEOD

Similar to the Indicators of the Peruvian Energy Sector, we first need to access the IEOD platform and check whether the download function uses POST or GET.

Figure 5: Figure 1

On the IEOD platform, select a date and the desired Excel file (e.g., “Anexo1_Resumen_0709.xlsx"). By right-clicking on the link, you can copy the URL address.

In this case, the copied address is: “https://www.coes.org.pe/portal/browser/download?url=Post%20Operaci%C3%B3n%2FReportes%2FIEOD%2F2020%2F09%20Setiembre%2F07%2FAnexo1_Resumen_0709.xlsx". From the bold portion, it is clear that this is a GET method (the data included in the request is visible in the URL, which is not the case with POST). However, what parts of the request change, and what do the strange characters like %2F, %20, or %C3%B3 mean?

By reviewing the Network section of the Chrome Developer Tools, we can see that the query starts after download?url=.

Figure 6: Figure 2

Based on the information above, we can infer that the GET portion carrying the data, i.e., the query, is:
*“Post Operación/Reportes/IEOD/2020/09 Setiembre/07/Anexo1_Resumen_**0709**.xlsx”*. Here’s what each bold element represents:

Year (e.g., 2020)
Month in numbers and letters (e.g., 09 Setiembre)
Day number (e.g., 07)
Day and month number combined (e.g., 0709)

If you click on view URL encoded, you will see that the encoded portion of the URL (Post%20Operaci%C3%B3n%2FReportes%2FIEOD%2F2020%2F09%20Setiembre%2F07%2FAnexo1_Resumen_0709.xlsx) translates as follows:

%20 = space
%C3%B3 = “ó” (accented letter)
%2F = “/” (slash)

In conclusion, the IEOD data is downloaded using the GET method, with the date details included in the request and embedded in the URL.

Downloading in R

With the knowledge gained earlier, we are ready to automate the download process. First, I will create the download URL by selecting the desired date and then converting it into the GET syntax identified earlier.

To work with dates, I will use the Lubridate library. To avoid potential issues with the spelling of the month “September” (Setiembre?), I will use a different date in October instead.

library(lubridate)

library(lubridate)
fecha = as.Date("2023/10/07")
Sys.setlocale("LC_TIME", "es_ES")

## [1] "es_ES"

format <- "%a@%A@%b@%B@%p@"
library(lubridate)
path = paste("https://www.coes.org.pe/portal/browser/download?url=Post%20Operaci%C3%B3n%2FReportes%2FIEOD%2F",
year(fecha), #El año
"%2F", # dash (/)
sprintf("%02d",month(fecha)), #número del mes, en 2 dígitos
"_", #espacio
str_to_title(month(fecha,label=TRUE,abbr = FALSE)), #mes en letras
#How to put the output of month() in spanish?
"%2F", # dash (/)
sprintf("%02d",day(fecha)), #número del día en 2 dígitos
"%2FAnexo1_Resumen_",
sprintf("%02d",day(fecha)), #numero del dia
sprintf("%02d",month(fecha)),#numero del mes
".xlsx",
sep="") #concateno las partes en uno solo
path

## [1] "https://www.coes.org.pe/portal/browser/download?url=Post%20Operaci%C3%B3n%2FReportes%2FIEOD%2F2023%2F10_Octubre%2F07%2FAnexo1_Resumen_0710.xlsx"

https://www.coes.org.pe/portal/browser/download?url=Post%20Operaci%C3%B3n%2FReportes%2FIEOD%2F2023%2F10_Octubre%2F19%2FAnexo1_Resumen_1910.xlsx Next, I will verify that the link works and successfully downloads electricity demand by geographic zone as well as demand for large free users (mining and manufacturing).

Downloading Geographic Zone Data

https://www.coes.org.pe/portal/browser/download?url=Post%20Operaci%C3%B3n%2FReportes%2FIEOD%2F2024%2F11_Noviembre%2F17%2FAnexoA_1711.xlsx

library(openxlsx)
zonas = read.xlsx(path, sheet = "DEMANDA_AREAS", rows = c(12, 64), cols = c(5:7)) / 1000
zonas

## NORTE CENTRO SUR
## 1 24.25947 94.15946 41.00531

As shown above, I downloaded the Excel file located at the path URL. From that Excel file, I selected the DEMANDA_AREAS sheet, which contains electricity demand in Peru, divided into three zones (North, Center, and South). From this sheet, I chose the columns (cols) containing demand by zone (columns 5 to 7). Additionally, I selected rows (rows) 12 (containing the column names) and 64 (containing the total daily consumption). Finally, I divided the values by 1,000 to convert the data from MWh to GWh.

Downloading Large Free Users Data (Mining and Manufacturing)

usuarios_libres = read.xlsx(path, sheet = "DEMANDA_UL", rows = c(20, 69))
duplicated_cols <- names(usuarios_libres)[duplicated(names(usuarios_libres))]
print(usuarios_libres)

## X1 CEMENTOS.PACASMAYO CREDITEX INDUSTRIAL.PIURA MINERA.HORIZONTE-LLACUABAMBA
## 1 MWh 1019.079 0 0 0
## YANACOCHA BARRICK.-.PIERINA BARRICK-.ALTO.CHICAMA MISKIMAYO-BAYOBAR
## 1 681.969 0 0 279.545
## GOLD.FIELD MARSA SIDER.PERU TRUPAL AGROPECUARIA.AURORA
## 1 629.1684 301.3809 979.8908 241.25 312.375
## PETROPERU.CL-6769.SE.GE1 PETROPERU.CL-6768.SE.GE1 ACEROS.AREQUIPA
## 1 587.759 0 1997.334
## JOCKEY.PLAZA ALICORP MINERA.ALPAYANA COMPAÑÍA.MINERA.CHUNGAR
## 1 0 229.4557 384.1315 337.6162
## COMPAÑÍA.MINERA.CONDESTABLE CREDITEX.-.INDEPENDENCIA CREDITEX.-.INGENIEROS
## 1 366.12 0 0
## EL.BROCAL ADMINISTRADORA.CERRO KIMBERLY.CLARK.PERU LECHE.GLORIA MARCOBRE
## 1 817.9798 93.22284 0 0 1153.913
## MEPSA MINERA.ANTAMINA MINERA.BUENAVENTURA.(UCHUCHACUA+HUANCAVELICA) CHINALCO
## 1 0 3036.311 122.29 3607.068
## MINERA.CORONA LOS.QUENUALES MINERA.SHOUXIN MINERA.VOLCAN MINSUR
## 1 131.204 0 324.9498 1141.029 67.78
## NEXA.RESOURCES.CAJAMARQUILLA NEXA.RESOURCES.-.DESIERTO NEXA.RESOURCES.-.PASCO
## 1 3910.832 773.12 328.9331
## OPP.FILM OXIDOS.PASCO PANASA SOFTYS.PERU QUIMPAC.-.PARAMONGA
## 1 76.87446 114.1832 104.4 278.3232 434.8822
## QUIMPAC.-.OQUENDO REFINERIA.LA.PAMPILLA SHOUGANG.HIERRO.PERU
## 1 1050.251 480.05 0
## TEJIDOS.SAN.JACINTO TRUPAL.-.ZARATE UNACEM-ANDINO UNACEM-LIMA MARCOBRE
## 1 109.318 209.0815 666.175 731.5305 1153.913
## ANGLO.AMERICAN.QUELLAVECO CAL.&.CEMENTO.SUR ANTAPACCAY HUDBAY.PERU
## 1 1254.671 74.97723 2740.34 0
## INDUSTRIAS.CACHIMAYO ARES-COTARUSE ARES-CALLALLI MINERA.BUENAVENTURA.-.ARES
## 1 442.9391 672.1465 18.83618 160.93
## MINERA.BUENAVENTURA.-.TAMBOMAYO CERRO.VERDE.(HIDROMETALÚRGICA)
## 1 176.635 610.3554
## CERRO.VERDE.-.SULFUROS CERRO.VERDE.(AMPLIACIÓN) LAS_BAMBAS
## 1 3371.444 6939.944 3678.651
## MINSUR.-.SAN.RAFAEL MINSUR.-.PUCAMARCA SPCC CEMENTOS.YURA
## 1 674.5289 0 6634.635 800.1874

# Make column names unique
names(usuarios_libres) <- make.names(names(usuarios_libres), unique = TRUE)
# Proceed with select() or other operations
usuarios_libres <- usuarios_libres %>% select(-X1)
sample(usuarios_libres, size = 4) # Display 4 sample companies

## SHOUGANG.HIERRO.PERU ALICORP MINERA.BUENAVENTURA...ARES HUDBAY.PERU
## 1 0 229.4557 160.93 0

As shown above, I downloaded the Excel file located at the path URL. From that Excel file, I selected the DEMANDA_UL sheet, which contains the demand data for large free users. From this sheet, I chose rows (rows) 20 (which contains company names) and 69 (which contains the total daily consumption).

Data Cleaning

Cleaning the name of the variables:

nombres = names(usuarios_libres) #nombre de las variables
nombres = gsub("\\."," ",nombres) #cambio puntos por espacios
nombres[which(nombres=="CERRO VERDE HIDROMETALÚRGICA")] = "REFINERIA CERRO VERDE"
nombres[which(nombres=="VOTORANTIM METAIS CAJAMARQUILLA")] = "REFINERIA CAJAMARQUILLA"
nombres[which(nombres=="ANTAMINA+09 OCTUBRE ")] = "ANTAMINA"
nombres[which(nombres=="UNACEM-LIMA")] = "CEMENTOS LIMA"
nombres[which(nombres=="UNACEM-ANDINO")] = "CEMENTOS ANDINO"
nombres[which(nombres=="INCASAC")] = "CACHIMAYO"
nombres[which(nombres=="SPCC")] = "SOUTHERN PERU"
names(usuarios_libres) = nombres
nombres

## [1] "CEMENTOS PACASMAYO"
## [2] "CREDITEX"
## [3] "INDUSTRIAL PIURA"
## [4] "MINERA HORIZONTE LLACUABAMBA"
## [5] "YANACOCHA"
## [6] "BARRICK PIERINA"
## [7] "BARRICK ALTO CHICAMA"
## [8] "MISKIMAYO BAYOBAR"
## [9] "GOLD FIELD"
## [10] "MARSA"
## [11] "SIDER PERU"
## [12] "TRUPAL"
## [13] "AGROPECUARIA AURORA"
## [14] "PETROPERU CL 6769 SE GE1"
## [15] "PETROPERU CL 6768 SE GE1"
## [16] "ACEROS AREQUIPA"
## [17] "JOCKEY PLAZA"
## [18] "ALICORP"
## [19] "MINERA ALPAYANA"
## [20] "COMPAÑÍA MINERA CHUNGAR"
## [21] "COMPAÑÍA MINERA CONDESTABLE"
## [22] "CREDITEX INDEPENDENCIA"
## [23] "CREDITEX INGENIEROS"
## [24] "EL BROCAL"
## [25] "ADMINISTRADORA CERRO"
## [26] "KIMBERLY CLARK PERU"
## [27] "LECHE GLORIA"
## [28] "MARCOBRE"
## [29] "MEPSA"
## [30] "MINERA ANTAMINA"
## [31] "MINERA BUENAVENTURA UCHUCHACUA HUANCAVELICA "
## [32] "CHINALCO"
## [33] "MINERA CORONA"
## [34] "LOS QUENUALES"
## [35] "MINERA SHOUXIN"
## [36] "MINERA VOLCAN"
## [37] "MINSUR"
## [38] "NEXA RESOURCES CAJAMARQUILLA"
## [39] "NEXA RESOURCES DESIERTO"
## [40] "NEXA RESOURCES PASCO"
## [41] "OPP FILM"
## [42] "OXIDOS PASCO"
## [43] "PANASA"
## [44] "SOFTYS PERU"
## [45] "QUIMPAC PARAMONGA"
## [46] "QUIMPAC OQUENDO"
## [47] "REFINERIA LA PAMPILLA"
## [48] "SHOUGANG HIERRO PERU"
## [49] "TEJIDOS SAN JACINTO"
## [50] "TRUPAL ZARATE"
## [51] "UNACEM ANDINO"
## [52] "UNACEM LIMA"
## [53] "MARCOBRE 1"
## [54] "ANGLO AMERICAN QUELLAVECO"
## [55] "CAL CEMENTO SUR"
## [56] "ANTAPACCAY"
## [57] "HUDBAY PERU"
## [58] "INDUSTRIAS CACHIMAYO"
## [59] "ARES COTARUSE"
## [60] "ARES CALLALLI"
## [61] "MINERA BUENAVENTURA ARES"
## [62] "MINERA BUENAVENTURA TAMBOMAYO"
## [63] "CERRO VERDE HIDROMETALÚRGICA "
## [64] "CERRO VERDE SULFUROS"
## [65] "CERRO VERDE AMPLIACIÓN "
## [66] "LAS_BAMBAS"
## [67] "MINSUR SAN RAFAEL"
## [68] "MINSUR PUCAMARCA"
## [69] "SOUTHERN PERU"
## [70] "CEMENTOS YURA"

The which function returns an index, which is the column number that satisfies the test (e.g., names == “X1”).

Some companies may have more than one measurement point. For instance, Cerro Verde has three: its old mine (Sulfuros), its mine expansion (Ampliación), and its hydrometallurgy center (Refinería). In this case, we sum the consumption of its old mine and the expansion, consolidating it into two measurement points (refinery and mine).

usuarios_libres = usuarios_libres %>%
mutate("CERRO VERDE" = `CERRO VERDE SULFUROS` + `CERRO VERDE AMPLIACIÓN `,
.after = `CERRO VERDE AMPLIACIÓN `) %>% # Sum two columns
select(-`CERRO VERDE SULFUROS`, -`CERRO VERDE AMPLIACIÓN `) # Remove old columns

With the data from the major electricity-demanding companies (large free users), I can categorize them based on their economic activities. The main sectors where electricity consumption is closely tied to production are mining and manufacturing.

I already know which columns correspond to mining companies, so I will select only those.

minerias = usuarios_libres %>%
select(1,2,4,5,6,9,11,12,17,20,24,25,27,28,30,31,32,33,34,35,38,40,41,42,43) # Select mining companies
names(minerias)

## [1] "CEMENTOS PACASMAYO"
## [2] "CREDITEX"
## [3] "MINERA HORIZONTE LLACUABAMBA"
## [4] "YANACOCHA"
## [5] "BARRICK PIERINA"
## [6] "GOLD FIELD"
## [7] "SIDER PERU"
## [8] "TRUPAL"
## [9] "JOCKEY PLAZA"
## [10] "COMPAÑÍA MINERA CHUNGAR"
## [11] "EL BROCAL"
## [12] "ADMINISTRADORA CERRO"
## [13] "LECHE GLORIA"
## [14] "MARCOBRE"
## [15] "MINERA ANTAMINA"
## [16] "MINERA BUENAVENTURA UCHUCHACUA HUANCAVELICA "
## [17] "CHINALCO"
## [18] "MINERA CORONA"
## [19] "LOS QUENUALES"
## [20] "MINERA SHOUXIN"
## [21] "NEXA RESOURCES CAJAMARQUILLA"
## [22] "NEXA RESOURCES PASCO"
## [23] "OPP FILM"
## [24] "OXIDOS PASCO"
## [25] "PANASA"

minerias %>%
transmute(Minas_GWh = rowSums(.) / 1000) # Sum mining company data

## Minas_GWh
## 1 17.71114

We can see that the mining companies consumed 2023-10-07 a total of 17.7 GWh.

Companies related to commerce/services, like Jockey Plaza Mall, do not have a strong relationship between production and electricity consumption, so they are excluded.

manufactureras = usuarios_libres %>%
select(-names(minerias)) # Select manufacturers
names(manufactureras)

## [1] "INDUSTRIAL PIURA" "BARRICK ALTO CHICAMA"
## [3] "MISKIMAYO BAYOBAR" "MARSA"
## [5] "AGROPECUARIA AURORA" "PETROPERU CL 6769 SE GE1"
## [7] "PETROPERU CL 6768 SE GE1" "ACEROS AREQUIPA"
## [9] "ALICORP" "MINERA ALPAYANA"
## [11] "COMPAÑÍA MINERA CONDESTABLE" "CREDITEX INDEPENDENCIA"
## [13] "CREDITEX INGENIEROS" "KIMBERLY CLARK PERU"
## [15] "MEPSA" "MINERA VOLCAN"
## [17] "MINSUR" "NEXA RESOURCES DESIERTO"
## [19] "SOFTYS PERU" "QUIMPAC PARAMONGA"
## [21] "QUIMPAC OQUENDO" "REFINERIA LA PAMPILLA"
## [23] "SHOUGANG HIERRO PERU" "TEJIDOS SAN JACINTO"
## [25] "TRUPAL ZARATE" "UNACEM ANDINO"
## [27] "UNACEM LIMA" "MARCOBRE 1"
## [29] "ANGLO AMERICAN QUELLAVECO" "CAL CEMENTO SUR"
## [31] "ANTAPACCAY" "HUDBAY PERU"
## [33] "INDUSTRIAS CACHIMAYO" "ARES COTARUSE"
## [35] "ARES CALLALLI" "MINERA BUENAVENTURA ARES"
## [37] "MINERA BUENAVENTURA TAMBOMAYO" "CERRO VERDE HIDROMETALÚRGICA "
## [39] "CERRO VERDE" "LAS_BAMBAS"
## [41] "MINSUR SAN RAFAEL" "MINSUR PUCAMARCA"
## [43] "SOUTHERN PERU" "CEMENTOS YURA"

manufactureras %>%
transmute(Manufact_GWh = rowSums(.) / 1000) # Calculate consumption

## Manufact_GWh
## 1 39.80477

We can see that the manufacturing companies consumed 2023-10-07 a total of 39.8 GWh.

Prefix	Symbol	`\(10^n\)`	Equivalent
kilo	K	`\(10^3\)`	1,000
mega	M	`\(10^6\)`	1,000,000
giga	G	`\(10^9\)`	1,000,000,000