BPA R Review

1. R Setup and Packages

Markdown Setup

The easiest way to manage packages is using pacman and its p_load() function.

Install pacman once in the RStudio Console. Then use p_load() to install missing packages and load them.

install.packages("pacman", repos = "https://cloud.r-project.org")

Package management functions

pacman::p_load()    # Installs missing packages and loads them
pacman::p_install() # Installs packages
pacman::p_unload()  # Unloads packages
pacman::p_update()  # Updates installed packages
pacman::p_version() # Checks package versions

Load packages

pacman::p_load(
  tidyverse,
  janitor,
  readxl,
  here,
  lubridate,
  skimr,
  broom,
  fixest
)

Package Reference

tidyverse — Data manipulation and visualization

A collection of packages for data manipulation, organization, visualization, and analysis.

Common functions:

select()       # Select columns
filter()       # Filter rows
mutate()       # Create or modify variables
summarise()    # Calculate summary statistics
group_by()     # Group observations
arrange()      # Sort observations
ggplot()       # Create visualizations

janitor — Data cleaning

Useful for cleaning messy datasets and creating frequency tables.

clean_names()       # Standardize column names
tabyl()             # Create frequency tables
get_dupes(data, ID) # Identify duplicate observations

readxl — Import Excel data

Used to import Excel spreadsheets into R.

read_excel("data.xlsx") # Import an Excel worksheet
excel_sheets("data.xlsx") # List worksheet names

here — Project file paths

Creates project-relative file paths rather than relying on absolute paths.

here("data", "raw", "data.csv")

lubridate — Dates and times

Simplifies working with dates, times, and time intervals.

year(date)              # Extract year
month(date)             # Extract month
day(date)               # Extract day
ymd("2026-11-12")       # Convert year-month-day string

skimr — Data summaries

Provides a detailed overview of a dataset, including descriptive statistics and missing values.

skim(data)

broom — Organize model results

Converts statistical model outputs into tidy data frames.

tidy(model)    # Coefficients, standard errors, p-values
glance(model)  # Model-level statistics
augment(model) # Observation-level results

fixest — Econometrics and regression

Useful for OLS, fixed-effects regression, instrumental variables, and models with robust or clustered standard errors.

feols(y ~ x1 + x2, data = data)

# Fixed effects
feols(y ~ x1 + x2 | year + region, data = data)

2. Importing and Inspecting Data

Importing and Exporting

read_csv("data.csv")       # Import CSV file
read_excel("data.xlsx")    # Import Excel file
readRDS("data.rds")        # Import R data object

write_csv(data, "data.csv") # Export CSV
saveRDS(data, "data.rds")   # Save R data object

Inspecting Data

head(data)       # Preview first observations
tail(data)       # Preview last observations
glimpse(data)    # Inspect columns and data types
dim(data)        # Number of rows and columns
nrow(data)       # Number of rows
ncol(data)       # Number of columns
names(data)      # Column names
summary(data)    # Descriptive statistics
skim(data)       # Detailed dataset overview
distinct(data)   # Identify unique observations

Common Data Types

Numeric

Numbers that can be used in mathematical calculations.

Examples: electricity price, load, temperature, revenue.

price <- 45.75
load_mw <- 1200

Integer

Whole numbers without decimal places.

year <- 2026L

Character

Text or strings.

region <- "Pacific Northwest"

Logical

Boolean values: TRUE or FALSE.

is_weekend <- TRUE

Date and datetime

Values representing calendar dates or specific times.

date <- as.Date("2026-11-12")

Factor

Categorical variables with defined levels.

region <- factor(c("North", "South", "North"))

Missing Values

Missing observations are represented by NA in R.

is.na(data)                 # Identify missing values
sum(is.na(data))            # Count missing values
colSums(is.na(data))        # Missing values by column
drop_na(data)               # Remove rows with missing values
replace_na(data, list(x=0)) # Replace missing values in x

Important: Missing values in calculations

Many R functions return NA when a vector contains missing values.

Use na.rm = TRUE to remove missing values from a calculation.

mean(data$price, na.rm = TRUE)
sum(data$load, na.rm = TRUE)

3. Data Manipulation

dplyr Functions

The dplyr package is part of the tidyverse and provides a consistent set of functions for manipulating datasets.

Select and filter — Choose columns and rows

select(data, price, load)       # Select columns
filter(data, region == "West")  # Select rows

• select() chooses which columns to keep.
• filter() chooses which observations to keep based on conditions.

Mutate and rename — Create and modify variables

mutate(data, revenue = price * load)

rename(data, electricity_price = price)

• mutate() creates new variables or modifies existing ones.
• rename() changes column names.

Arrange — Sort observations

arrange(data, date)          # Ascending
arrange(data, desc(price))   # Descending

Group and summarize — Aggregate observations

data |>
  group_by(region) |>
  summarise(
    average_price = mean(price, na.rm = TRUE),
    total_load = sum(load, na.rm = TRUE),
    observations = n()
  )

• group_by() defines groups for subsequent calculations.
• summarise() calculates summary statistics.
• n() counts observations within each group.
• count() is a convenient shortcut for counting observations by group.

Other useful dplyr functions

count(data, region)          # Count observations
distinct(data, region)       # Unique values or combinations
across()                     # Apply functions across columns
case_when()                  # Create conditional variables
lag(price)                   # Previous observation
lead(price)                  # Next observation

Combining and Reshaping Datasets

Joining datasets

Joins combine datasets using one or more shared identifiers.

left_join(x, y, by = "ID")   # Keep all rows in x
right_join(x, y, by = "ID")  # Keep all rows in y
inner_join(x, y, by = "ID")  # Keep matching rows only
full_join(x, y, by = "ID")   # Keep all rows from both
anti_join(x, y, by = "ID")   # Keep unmatched rows from x

Important: Always check whether join keys are unique before joining. Duplicate keys can unexpectedly multiply rows and affect subsequent statistics and regressions.

count(data, ID)
janitor::get_dupes(data, ID)

Reshaping datasets

Reshaping changes how observations and variables are organized.

pivot_longer(
  data,
  cols = c(January, February, March),
  names_to = "month",
  values_to = "sales"
)

pivot_wider(
  data,
  names_from = month,
  values_from = sales
)

• pivot_longer() converts wide data into long data by stacking columns into key-value rows.
• pivot_wider() converts long data into wide data by spreading values across columns.

Pipes

Pipes pass the result of one operation into the next operation.

data |>
  filter(region == "West") |>
  group_by(year) |>
  summarise(
    average_price = mean(price, na.rm = TRUE)
  )

Common pipe operators:

• |> Base R pipe.
• %>% magrittr pipe, commonly used in tidyverse code.

Both allow operations to be written sequentially rather than nesting functions.

4. Statistics and Econometrics

Descriptive Statistics

Descriptive Statistics

Numerical summaries that describe the characteristics, distribution, and variability of a dataset.

Measure	Definition
Mean	Arithmetic average
Median	Middle observation
Variance	Average squared deviation from the mean
Standard deviation	Square root of variance; measures spread in original units
Quantiles	Values dividing a distribution into specified proportions
Interquartile range (IQR)	75th percentile minus 25th percentile
Minimum / Maximum	Smallest and largest observations

Why they matter: Help identify unusual values, skewed distributions, missing data, and differences between groups.

R functions:

summary(data)
mean(x, na.rm = TRUE)
median(x, na.rm = TRUE)
sd(x, na.rm = TRUE)
var(x, na.rm = TRUE)
quantile(x, na.rm = TRUE)

Correlations

Measures the strength and direction of the relationship between two variables.

Pearson’s correlation coefficient, r, ranges from -1 to +1.

• +1: Perfect positive linear relationship.
• 0: No linear relationship.
• -1: Perfect negative linear relationship.

Why they matter: Useful for exploring relationships between variables and identifying potential explanatory variables before regression analysis.

Important: Correlation does not imply causation. A correlation may reflect reverse causality, a third variable, or a shared trend.

R functions:

cor(x, y)
cor.test(x, y)
GGally::ggpairs(data)

Linear Regression

A statistical method used to estimate the relationship between an outcome variable (dependent variable) and one or more explanatory variables (independent variables).

Simple linear regression:

\[ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i \]

Where:

• Y: Outcome variable.
• X: Explanatory variable.
• β₀: Intercept.
• β₁: Slope coefficient.
• ε: Error term.

Interpretation: A coefficient represents the estimated change in the outcome associated with a one-unit increase in an explanatory variable, holding other included variables constant.

Example:

If the estimated coefficient on electricity load is 0.04, a one-unit increase in load is associated with a 0.04-unit increase in price, holding other included variables constant.

R functions:

lm(y ~ x, data = data)

lm(y ~ x1 + x2 + x3, data = data)

fixest::feols(y ~ x1 + x2, data = data)

Model Summaries

Regression output summarizes coefficient estimates, uncertainty, and model fit.

Output	Meaning
Coefficient	Estimated relationship between an explanatory variable and outcome
Standard error	Estimated uncertainty around a coefficient
t-statistic	Coefficient divided by its standard error
p-value	Evidence against the null hypothesis under the model assumptions
R-squared	Proportion of sample variation in the outcome explained by the fitted model
Adjusted R-squared	R-squared adjusted for the number of explanatory variables
F-statistic	Tests a joint restriction on model coefficients
Observations	Number of observations used in estimation

Important: A high R-squared does not automatically mean a model is unbiased, causal, or useful for forecasting.

R functions:

summary(model)
broom::tidy(model)
broom::glance(model)

Predictions

Estimated outcome values generated using a fitted statistical model and specified explanatory-variable values.

Two important distinctions:

• Fitted values: Predictions for observations already included in the estimation dataset.
• Out-of-sample predictions: Predictions for new observations not used to estimate the model.

Why they matter: Useful for forecasting, scenario analysis, and estimating expected outcomes under different conditions.

Important: A prediction is not necessarily a causal estimate.

R functions:

predict(model)
predict(model, newdata = new_data)

Residuals

The differences between observed outcomes and the values predicted by a regression model.

\[ e_i = Y_i - \hat{Y}_i \]

• Positive residual: Model underpredicted the observed value.
• Negative residual: Model overpredicted the observed value.
• Zero residual: Prediction exactly matches the observation.

Why they matter: Residuals help evaluate model fit and identify potential problems such as nonlinear relationships, outliers, heteroskedasticity, and autocorrelation.

R functions:

residuals(model)
fitted(model)
broom::augment(model)
plot(model)

Hypothesis Testing

A statistical procedure used to evaluate whether sample data provide sufficient evidence against a specified claim about a population or model parameter.

Two competing hypotheses:

• Null hypothesis (H₀): The default claim, often that a coefficient equals zero or that no relationship exists.
• Alternative hypothesis (H₁): The competing claim, such as a coefficient being different from zero.

Important terminology:

• p-value: Probability of observing a test statistic at least as extreme as the one obtained, assuming the null hypothesis and model assumptions are true.
• Significance level (α): Threshold chosen for rejecting the null hypothesis, commonly 0.05.
• Type I error: Rejecting a true null hypothesis.
• Type II error: Failing to reject a false null hypothesis.

Important: Statistical significance does not necessarily imply practical importance. We generally say “reject” or “fail to reject” the null hypothesis rather than prove or accept it.

R functions:

t.test(x, y)
summary(model)

Robust Standard Errors

Standard errors calculated using an estimator that is robust to certain violations of conventional regression assumptions, particularly heteroskedasticity.

Why they matter: Conventional standard errors may be incorrect when error variance is not constant, leading to unreliable t-statistics, p-values, and confidence intervals.

Important distinctions:

• Regression coefficients generally remain the same.
• Standard errors change.
• t-statistics, p-values, and confidence intervals may change.

Robust standard errors do not fix biased coefficients, omitted variables, or endogeneity.

R functions:

sandwich::vcovHC(model, type = "HC1")

lmtest::coeftest(
  model,
  vcov = sandwich::vcovHC(model, type = "HC1")
)

Clustered Standard Errors

Standard errors that account for the possibility that observations within the same group are correlated.

Example: Hourly electricity prices within a region may share common influences such as weather, transmission infrastructure, and generation resources.

Clustering by region allows errors to be correlated within each region while assuming independence across regions.

Why they matter: Treating correlated observations as independent can underestimate uncertainty and produce misleading statistical significance.

Important: The number of independent clusters matters. A large number of observations does not compensate for having very few clusters.

R functions:

fixest::feols(
  price ~ load + temperature,
  data = energy_data,
  vcov = ~region
)

Fixed Effects and Panel Regression

Fixed effects control for unobserved characteristics of individual units that do not change over time.

Panel data follow the same units (such as regions, states, or firms) across multiple time periods.

A simplified model:

\[ Y_{it} = \beta X_{it} + \alpha_i + \epsilon_{it} \]

Where:

• Y: Outcome for unit i at time t.
• X: Explanatory variable.
• β: Estimated coefficient.
• α: Unit-specific fixed effect.
• ε: Error term.

Interpretation: The coefficient is identified using variation within the same unit over time, rather than simply comparing differences between units.

Important: Time-invariant explanatory variables are absorbed by unit fixed effects. Fixed effects do not automatically control for unobserved factors that change over time.

R functions:

fixest::feols(
  price ~ load + temperature | region,
  data = energy_data
)

# Region and year fixed effects
fixest::feols(
  price ~ load + temperature | region + year,
  data = energy_data
)

Weighted Regression

A regression technique that assigns different weights to observations when estimating the model.

Ordinary least squares gives every observation equal weight. Weighted least squares gives observations different influence on the estimated coefficients.

\[ \min_{\beta} \sum_i w_i(Y_i-X_i\beta)^2 \]

Weights may be used for:

• Unequal sampling probabilities.
• Different levels of measurement precision.
• Appropriately modeled heteroskedasticity.

Important: Survey weights, frequency weights, and inverse-variance weights have different meanings and should not be used interchangeably.

R functions:

lm(
  price ~ load + temperature,
  data = energy_data,
  weights = precision_weight
)

Instrumental Variables (IV)

A method used to estimate causal relationships when an explanatory variable is endogenous.

Endogeneity occurs when an explanatory variable is correlated with the regression error term.

Potential causes include:

• Reverse causality.
• Omitted variable bias.
• Measurement error.
• Simultaneous determination of variables.

An instrument must satisfy two key conditions:

1. Relevance: The instrument is correlated with the endogenous explanatory variable.
2. Exogeneity / exclusion restriction: The instrument affects the outcome only through the endogenous explanatory variable and is uncorrelated with the structural error term.

Example: Electricity prices and consumption may be simultaneously determined. An external source of price variation might help identify the effect of prices on consumption, provided the instrument is valid.

Important: A strong first-stage relationship does not prove an instrument is valid. Weak instruments can produce unreliable estimates, and the exclusion restriction is generally not directly testable.

R functions:

fixest::feols(
  consumption ~ temperature | price ~ instrument,
  data = energy_data
)

Heteroskedasticity

A condition in which the variance of the regression error term is not constant across observations.

• Homoskedasticity: Constant error variance.
• Heteroskedasticity: Nonconstant error variance.

Why it matters: OLS coefficients can remain unbiased under the usual exogeneity assumptions, but conventional standard errors are generally incorrect.

How to address it:

• Use heteroskedasticity-robust standard errors.
• Consider weighted least squares when the variance structure is appropriately understood.
• Reconsider model specification if necessary.

R functions:

lmtest::bptest(model)
sandwich::vcovHC(model, type = "HC1")

Multicollinearity

A condition in which two or more explanatory variables are highly correlated with one another.

Why it matters:

• Coefficient estimates may have large standard errors.
• Individual coefficients may appear statistically insignificant.
• Coefficient estimates may be sensitive to model specification.
• It becomes difficult to isolate individual relationships.

Important: Multicollinearity does not automatically bias OLS estimates. It primarily affects precision and stability.

How to identify it:

• Examine correlation matrices.
• Calculate variance inflation factors (VIFs).

R functions:

cor(data)
car::vif(model)

Omitted Variable Bias (OVB)

Bias in regression coefficient estimates that occurs when an important variable is excluded from a model and is correlated with both an included explanatory variable and the outcome.

Two conditions:

1. The omitted variable affects the outcome.
2. The omitted variable is correlated with an included explanatory variable.

Example: Estimating the relationship between electricity prices and consumption without accounting for seasonality may bias the estimated price relationship if seasonality affects both.

How to address it:

• Include relevant control variables.
• Use fixed effects when appropriate.
• Consider instrumental variables when endogeneity is present and a credible instrument exists.
• Improve the research design.

Important: Adding more variables does not automatically improve a model. Inappropriate controls can introduce other forms of bias.

Fixed Effects vs. Random Effects

Two approaches to modeling differences between individual units in panel data.

Feature	Fixed Effects	Random Effects
Unit-specific differences	May be correlated with explanatory variables	Assumed uncorrelated with explanatory variables
Estimation	Uses within-unit variation	Uses within- and between-unit variation
Time-invariant variables	Cannot estimate separate coefficients	Can estimate coefficients
Main assumption	Allows correlation between unit effects and regressors	Assumes no correlation

Interpretation:

• Fixed effects control for unobserved unit characteristics that may be correlated with explanatory variables.
• Random effects assume those characteristics are uncorrelated with explanatory variables.

How to choose: Fixed effects are often preferable when unobserved unit-level characteristics may be correlated with regressors. Random effects may be appropriate when the stronger assumption is defensible.

R functions:

plm::plm(
  y ~ x1 + x2,
  data = panel_data,
  index = c("region", "year"),
  model = "within"
)

plm::plm(
  y ~ x1 + x2,
  data = panel_data,
  index = c("region", "year"),
  model = "random"
)

plm::phtest(fixed_model, random_model)

Interaction Terms

Variables created by multiplying two or more explanatory variables together. Interactions allow the relationship between one variable and the outcome to depend on another variable.

Example: The relationship between temperature and electricity demand may differ between weekdays and weekends.

\[ Load_i = \beta_0 + \beta_1 Temperature_i + \beta_2 Weekend_i + \beta_3(Temperature_i \times Weekend_i) + \epsilon_i \]

Interpretation:

• β₀: Expected load when temperature is zero on weekdays.
• β₁: Temperature slope on weekdays.
• β₂: Weekend versus weekday difference when temperature is zero.
• β₃: Difference in the temperature slope between weekends and weekdays.

The temperature slope is:

• Weekdays: β₁
• Weekends: β₁ + β₃

R syntax:

lm(
  load ~ temperature * weekend,
  data = energy_data
)

The * operator includes both main effects and the interaction:

temperature * weekend

# Equivalent to:
temperature + weekend + temperature:weekend

5. Data Visualization

ggplot2

The ggplot2 package creates visualizations using a layered approach.

Basic structure:

ggplot(data, aes(x = variable1, y = variable2)) +
  geom_point() +
  labs(
    title = "Relationship Between Variables",
    x = "Variable 1",
    y = "Variable 2"
  )

Common Geometries

Scatterplots

Show the relationship between two continuous variables.

ggplot(data, aes(x = load, y = price)) +
  geom_point()

Line graphs

Useful for visualizing trends over time.

ggplot(data, aes(x = date, y = price)) +
  geom_line()

Bar charts

Compare categories or groups.

ggplot(data, aes(x = region, y = average_price)) +
  geom_col()

Histograms

Visualize the distribution of a continuous variable.

ggplot(data, aes(x = price)) +
  geom_histogram()

6. Time-Series Analysis

Dates and Time

Electricity-market data frequently involve hourly, daily, monthly, and annual observations.

ymd("2026-11-12")       # Parse a date
year(date)              # Extract year
month(date)             # Extract month
day(date)               # Extract day
hour(datetime)          # Extract hour
floor_date(date, "month") # Round down to month

Time-Series Concepts

Trends

Long-term increases or decreases in a variable over time.

Seasonality

Recurring patterns that occur at regular intervals, such as daily, weekly, or annual cycles.

Autocorrelation

The correlation of a variable with its own past values.

Lagged Variables

Variables representing previous observations.

data |>
  arrange(date) |>
  mutate(previous_price = lag(price))

Moving Averages

Averages calculated over a rolling window of observations to smooth short-term fluctuations.

7. Electricity Markets and Economic Applications

Electricity Market Fundamentals

Load and Demand

Electricity consumption over a specified period, commonly measured in megawatts (MW).

Generation Mix

The combination of resources used to generate electricity, including hydropower, natural gas, coal, wind, solar, and nuclear.

Marginal Cost

The additional cost of producing one more unit of electricity.

Locational Marginal Pricing (LMP)

The price of electricity at a particular location, reflecting the marginal cost of serving additional demand, including energy, congestion, and losses.

Transmission Congestion

Occurs when transmission constraints limit the ability to move electricity between locations.

Resource Adequacy

The ability of the electricity system to maintain sufficient resources to meet demand reliably.

Ancillary Services

Services that help maintain grid reliability, including frequency regulation and operating reserves.

Hydropower

Electricity generated from moving water. Reservoir levels, streamflow, water availability, and operating constraints influence hydropower generation.

Western Electricity Markets

CAISO

California Independent System Operator.

WEIM

Western Energy Imbalance Market.

EDAM

Extended Day-Ahead Market.

Markets+

A Western day-ahead electricity market initiative.

RTO / ISO

Regional Transmission Organization / Independent System Operator. Organizations that coordinate electricity transmission and wholesale market operations.

8. Useful R Shortcuts and Troubleshooting

Common RStudio Shortcuts

Action	Mac	Windows / Linux
Run current line or selection	Cmd + Enter	Ctrl + Enter
Insert assignment operator	Option + -	Alt + -
Reformat code	Cmd + Shift + A	Ctrl + Shift + A
Comment / uncomment	Cmd + Shift + C	Ctrl + Shift + C
Insert code chunk	Option + Cmd + I	Ctrl + Alt + I

Common Troubleshooting

Object not found

R cannot find the object referenced in the code.

Check spelling, capitalization, and whether the object has been created or imported.

Could not find function

The function may belong to a package that has not been loaded.

Load the package or use its namespace:

library(dplyr)

dplyr::filter(data, price > 50)

Unexpected NA values

Check whether missing observations or mismatched data types are affecting calculations and joins.

sum(is.na(data))
str(data)

Unexpected number of rows after joining

Check whether the join keys contain duplicates.

count(data, ID)
janitor::get_dupes(data, ID)

Package installation errors

Make sure R has a CRAN repository selected.

install.packages(
  "package_name",
  repos = "https://cloud.r-project.org"
)

9. Experimental Data Analysis

Importing and defining data

Data downloaded from https://transmission.bpa.gov/business/operations/wind/baltwg3.aspx. Real-time data from the last seven days, updated every five minutes.

Includes data on load and VER (Variable Energy Resource, usually solar or wind, which change with the conditions), as well as hydro, fossil/biomass, interchange, and nuclear.

Term	Definition	Constraints
Load	The amount of electrical power being consumed by customers and other electricity users at any given point in time. Usually measured in megawatts (MW). This represents electricity demand not generation.
VER (Variable Energy Resources)	Electricity generation resources whose output varies based on environmental conditions rather than being fully dispatchable on demand. Common examples include wind or solar generation.	Dependent on weather conditions, which can be unstable, but not always unpredictable.
Hydroelectric	Generated by converting the energy of moving or falling water into electrical energy using turbines and generators. Can be adjusted relatively quickly by controlling water flowing through turbines.	Water availability, reservoir levels, environmental requirements, and operational restrictions.
Fossil/Biomass Generation	Generated by the burning of fossil fuels of biomass. Examples: Natural gas, coal, organic matter used as fuel. Most fossil fuel plants are dispatchable, and output is easily adjusted.	Emissions and environmental regulations. Operating costs depend on fuel prices.
Nuclear	Generated through nuclear fission, where splitting atomic nuclei releases heat which produces steam to drive a turbine. Provides relatively consistent, continuous generation, but has high construction cost and is generally less flexible.
Net Interchange	The net amount of electricity transferred between a balancing authority (ie BPA) and neighboring balancing authorities.
Generation	Electricity produced by power plants and other resources.
Dispatch	The process of determining which generation resources operate and at what output levels.
Dispatchable Generation	Generation that operators can control or adjust when needed, subject to operating constraints.
Ramp Rate	How quickly a generation resource can increase or decrease its output.
Peak Load	The maximum electricity demand during a specific period.
Load Forecasting	Predicting the future electricity demand.
Generation Mix	The combination of different resources used to produce electricity.
MW (megawatt)	A unit of power, measuring the instantaneous rate of electricity production or consumption.
MWh (megawatt-hour)	A unit of energy, measuring electricity produced or consumed over time.

Note: The dataset reports five-minute MW readings. These are measurements of power at particular points in time, not the total energy generated over each five-minute interval.

Simply, \[Generation + Net Imports ≈ Load \] Expanding on that, \[Hydro + VER + Fossil/Biomass + Nuclear + Net Imports ≈ Load \]

# Import BPA balancing authority data
load_data <- read_tsv(
  "~/Desktop/BPA_Data.txt"
) %>%
  clean_names() %>%
  mutate(
    date_time = mdy_hm(date_time, tz = "America/Los_Angeles")
  )

# Inspect imported data
# glimpse(load_data)
# summary(load_data)

Exploratory data analysis

Data overview & cleaning

# Number of observations and variables
dim(load_data)

#> [1] 2016    7

# Variable names
names(load_data)

#> [1] "date_time"      "load"           "ver"            "hydro"         
#> [5] "fossil_biomass" "interchange"    "nuclear"

# Check data types
glimpse(load_data)

#> Rows: 2,016
#> Columns: 7
#> $ date_time      <dttm> 2026-11-07 00:00:00, 2026-11-07 00:05:00, 2026-11-07 0…
#> $ load           <dbl> 6191, 6188, 6189, 6172, 6142, 6143, 6133, 6095, 6030, 6…
#> $ ver            <dbl> 446, 391, 373, 361, 363, 363, 372, 384, 380, 379, 380, …
#> $ hydro          <dbl> 5484, 5547, 5502, 5511, 5517, 5484, 5467, 5508, 5458, 5…
#> $ fossil_biomass <dbl> 779, 781, 781, 782, 782, 782, 782, 780, 781, 782, 781, …
#> $ interchange    <dbl> 1702, 1717, 1656, 1674, 1699, 1671, 1677, 1766, 1774, 1…
#> $ nuclear        <dbl> 1144, 1146, 1148, 1151, 1149, 1146, 1149, 1151, 1145, 1…

# Check for missing values
colSums(is.na(load_data))

#>      date_time           load            ver          hydro fossil_biomass 
#>              0            159            159            159            159 
#>    interchange        nuclear 
#>            159            159

# Check for duplicate timestamps
sum(duplicated(load_data$date_time))

#> [1] 0

# Check date range
range(load_data$date_time)

#> [1] "2026-11-07 00:00:00 PST" "2026-11-13 23:55:00 PST"

# We need to drop for November 13, since we only have partial day data. 

# Count missing observations by date
load_data %>%
  group_by(date = as.Date(date_time)) %>%
  summarise(
    total_observations = n(),
    missing_load = sum(is.na(load)),
    valid_load = sum(!is.na(load)),
    .groups = "drop"
  )

# Display timestamps with missing measurements
load_data %>%
  filter(is.na(load)) %>%
  select(date_time)

  # print(n = 30)

# Retain only complete calendar days
complete_data <- load_data %>%
  mutate(date = as.Date(date_time)) %>%
  group_by(date) %>%
  filter(
    n() == 288,
    all(
      if_all(
        c(load, ver, hydro, fossil_biomass, interchange, nuclear),
        ~ !is.na(.x)
      )
    )
  ) %>%
  ungroup()

# Calculate descriptive statistics
descriptive_stats <- complete_data %>%
  summarise(
    across(
      c(load, ver, hydro, fossil_biomass, interchange, nuclear),
      list(
        N = ~sum(!is.na(.x)),
        Mean = ~mean(.x, na.rm = TRUE),
        Median = ~median(.x, na.rm = TRUE),
        SD = ~sd(.x, na.rm = TRUE),
        Min = ~min(.x, na.rm = TRUE),
        Max = ~max(.x, na.rm = TRUE)
      )
    )
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = c("Variable", ".value"),
    names_pattern = "(.+)_(N|Mean|Median|SD|Min|Max)"
  ) %>%
  mutate(
    Variable = recode(
      Variable,
      load = "Load",
      ver = "Variable Energy Resources",
      hydro = "Hydro",
      fossil_biomass = "Fossil/Biomass",
      interchange = "Net Interchange",
      nuclear = "Nuclear"
    ),
    across(where(is.numeric), ~round(.x, 2))
  )

# Display formatted table
descriptive_stats %>%
  knitr::kable(
    caption = "Descriptive Statistics: BPA Balancing Authority",
    col.names = c(
      "Variable",
      "N",
      "Mean (MW)",
      "Median (MW)",
      "Std. Dev. (MW)",
      "Minimum (MW)",
      "Maximum (MW)"
    ),
    align = c("l", rep("r", 6))
  ) %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE,
    position = "center"
  )

Descriptive Statistics: BPA Balancing Authority

Variable	N	Mean (MW)	Median (MW)	Std. Dev. (MW)	Minimum (MW)	Maximum (MW)
Load	1440	6896.68	6955.5	590.26	5732	8132
Variable Energy Resources	1440	414.26	161.0	534.41	5	2367
Hydro	1440	6567.88	6539.0	968.35	4058	8789
Fossil/Biomass	1440	984.71	1092.0	253.63	489	1200
Net Interchange	1440	2258.71	2211.5	518.70	1082	3790
Nuclear	1440	1147.98	1148.0	4.00	1125	1157

# Summarize data completeness by day
daily_data_quality <- load_data %>%
  mutate(date = as.Date(date_time)) %>%
  group_by(date) %>%
  summarise(
    observations = n(),
    missing = sum(is.na(load)),
    valid = sum(!is.na(load)),
    completeness = round(valid / 288 * 100, 1),
    .groups = "drop"
  )

# Display data-quality table
daily_data_quality %>%
  knitr::kable(
    caption = "Daily Data Completeness",
    col.names = c(
      "Date",
      "Observations",
      "Missing",
      "Valid",
      "Completeness (%)"
    ),
    align = c("l", "r", "r", "r", "r")
  ) %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE,
    position = "center"
  )

Daily Data Completeness

Date	Observations	Missing	Valid	Completeness (%)
2026-11-07	192	0	192	66.7
2026-11-08	288	0	288	100.0
2026-11-09	288	0	288	100.0
2026-11-10	288	0	288	100.0
2026-11-11	288	0	288	100.0
2026-11-12	288	0	288	100.0
2026-11-13	288	63	225	78.1
2026-11-14	96	96	0	0.0

Daily summaries: load and generation

For this, i will ise complete_data so we are comparing full days rather than accidentally including incomplete days. For each day, we will calculate:

• Average Load: typical balancing area load throughout the day
• Peak Load: highest five-minute load reading
• Minimum load: lowest five-minute load reading
• Average generation by resource
• Average interchange

# Calculate daily load and generation summaries
daily_summary <- complete_data %>%
  group_by(date) %>%
  summarise(
    # Date
    day_of_week = wday(first(date), label = TRUE, abbr = FALSE), #extracts day of week, returning teh first value in the "date" group

    # Load
    observations = n(),
    avg_load = mean(load),
    peak_load = max(load),
    peak_time = format(date_time[which.max(load)], "%H:%M"),
    min_load = min(load),
    min_time = format(date_time[which.min(load)], "%H:%M"),

    # Generation
    avg_ver = mean(ver),
    avg_hydro = mean(hydro),
    avg_fossil_biomass = mean(fossil_biomass),
    avg_nuclear = mean(nuclear),

    # Interchange
    avg_interchange = mean(interchange),

    .groups = "drop"
  )

# Display daily summary table
daily_summary %>%
  knitr::kable(
    caption = "Daily Load and Generation Summary",
    col.names = c(
      "Date",
      "Day",
      "Observations",
      "Average Load (MW)",
      "Peak Load (MW)",
      "Peak Time",
      "Minimum Load (MW)",
      "Minimum Time",
      "Average VER (MW)",
      "Average Hydro (MW)",
      "Average Fossil/Biomass (MW)",
      "Average Nuclear (MW)",
      "Average Interchange (MW)"
    ),
    align = c("l", "l", rep("r", 11)),
    digits = 1
  ) %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE,
    position = "center"
  )

Daily Load and Generation Summary

Date	Day	Observations	Average Load (MW)	Peak Load (MW)	Peak Time	Minimum Load (MW)	Minimum Time	Average VER (MW)	Average Hydro (MW)	Average Fossil/Biomass (MW)	Average Nuclear (MW)	Average Interchange (MW)
2026-11-08	Sunday	288	6481.1	7057	17:35	5732	02:00	1240.4	5813.8	549.7	1144.6	2308.1
2026-11-09	Monday	288	7009.9	7779	07:40	5979	02:00	125.6	6790.7	985.6	1150.4	2082.6
2026-11-10	Tuesday	288	7206.6	8132	17:30	6030	03:00	523.8	6874.0	1142.9	1149.4	2524.0
2026-11-11	Wednesday	288	6872.9	7820	17:30	5931	03:30	128.3	6667.7	1112.2	1147.5	2223.8
2026-11-12	Thursday	288	6912.9	7581	08:10	5893	02:30	53.1	6693.2	1133.1	1148.0	2155.0

The above table shows all of the complete days in the dataset (November 8 - November 12). The lowest load occurred on November 8, which is s Sunday. This tracks, as industry is not as active on weekends. We see a higher interchange on this day as well. The minimum load time is relatively consistent across each day, ranging from 2am - 3:30am. Again, this makes sense, as the vast majority of people are asleep at this time.

We see the peak load times vary, however. Two days (the 9th and the 12th) see peak times in the morning, at 7:40am and 8:10am respectively. The other days have almost remarkable consistency, with peak load occurring in the five minutes between 5:30pm and 5:35pm.

The average VER is significantly higher on Sunday. Sunday was sunnier than other days, so this may indicate a higher access to solar power.

Daily energy calculations: MW to MWh

Remember: MW (megawatts) are a measure of power at a point in time. MWh (megawatt-hours) are a measure of energy over a period of time. The data contains five-minute MW observations.

We can estimate energy by multiplying each observation by the length of its interval: \[ MWh = MW \times \frac{5}{60}\]

We can now calculate daily energy using complete_data again.

# Calculate daily energy consumption and generation
daily_energy <- complete_data %>%
  group_by(date) %>%
  summarise(
    # Energy consumption
    load_mwh = sum(load * (5 / 60)),

    # Energy generation by resource
    ver_mwh = sum(ver * (5 / 60)),
    hydro_mwh = sum(hydro * (5 / 60)),
    fossil_biomass_mwh = sum(fossil_biomass * (5 / 60)),
    nuclear_mwh = sum(nuclear * (5 / 60)),

    # Total reported generation
    total_generation_mwh = sum(
      (ver + hydro + fossil_biomass + nuclear) * (5 / 60)
    ),

    .groups = "drop"
  )

# Display daily energy table
daily_energy %>%
  knitr::kable(
    caption = "Daily Energy Consumption and Generation",
    col.names = c(
      "Date",
      "Load (MWh)",
      "VER (MWh)",
      "Hydro (MWh)",
      "Fossil/Biomass (MWh)",
      "Nuclear (MWh)",
      "Total Generation (MWh)"
    ),
    align = c("l", rep("r", 6)),
    digits = 1
  ) %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE,
    position = "center"
  )

Daily Energy Consumption and Generation

Date	Load (MWh)	VER (MWh)	Hydro (MWh)	Fossil/Biomass (MWh)	Nuclear (MWh)	Total Generation (MWh)
2026-11-08	155545.8	29770.3	139530.7	13192.5	27470.9	209964.4
2026-11-09	168237.1	3013.7	162977.8	23655.3	27608.4	217255.2
2026-11-10	172957.7	12571.7	164976.2	27429.9	27584.8	232562.8
2026-11-11	164950.7	3079.4	160024.7	26692.8	27541.0	217338.0
2026-11-12	165910.8	1275.6	160636.2	27195.0	27552.1	216658.9

Hourly load profiles

First we need to calculate the average load by hour. Each hour contains 12 5-minute observations. We will group all five days together by hour, then calculate the average load across all observations recorded during that hour.

# Calculate average load by hour of day
hourly_load <- complete_data %>%
  mutate(
    hour = lubridate::hour(date_time)
  ) %>%
  group_by(hour) %>%
  summarise(
    avg_load = mean(load),
    min_load = min(load),
    max_load = max(load),
    sd_load = sd(load),
    observations = n(),
    .groups = "drop"
  )

# Display hourly load profile
hourly_load %>%
  knitr::kable(
    caption = "Average Hourly BPA Balancing Authority Load",
    col.names = c(
      "Hour",
      "Average Load (MW)",
      "Minimum Load (MW)",
      "Maximum Load (MW)",
      "Std. Dev. (MW)",
      "Observations"
    ),
    align = c("l", rep("r", 5)),
    digits = 1
  ) %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE,
    position = "center"
  )

Average Hourly BPA Balancing Authority Load

Hour	Average Load (MW)	Minimum Load (MW)	Maximum Load (MW)	Std. Dev. (MW)	Observations
0	6129.3	5837	6421	138.9	60
1	6021.7	5738	6271	134.8	60
2	5978.7	5732	6184	122.0	60
3	5997.3	5742	6178	127.1	60
4	6162.4	5790	6503	197.9	60
5	6528.3	5894	7030	318.6	60
6	7005.9	6142	7629	428.2	60
7	7306.6	6396	7779	427.4	60
8	7364.9	6728	7778	319.7	60
9	7354.1	6869	7639	265.2	60
10	7291.5	6841	7641	257.0	60
11	7232.4	6822	7599	264.3	60
12	7199.9	6759	7575	283.2	60
13	7161.6	6701	7548	296.1	60
14	7116.7	6659	7485	300.5	60
15	7154.4	6678	7655	308.4	60
16	7268.9	6606	8030	396.0	60
17	7526.4	6973	8132	362.9	60
18	7462.7	6875	8062	353.6	60
19	7334.5	6773	7949	344.0	60
20	7150.4	6600	7734	300.0	60
21	6885.8	6415	7434	256.8	60
22	6578.3	6180	7078	208.9	60
23	6307.5	5998	6637	161.0	60

# Plot average hourly load
ggplot(hourly_load, aes(x = hour, y = avg_load)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_x_continuous(
    breaks = 0:23,
    labels = sprintf("%02d:00", 0:23)
  ) +
  labs(
    title = "Average Hourly BPA Balancing Authority Load",
    subtitle = "Average five-minute load observations across November 8–12, 2026",
    x = "Hour of Day",
    y = "Average Load (MW)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(face = "bold")
  )

This is the pooled average, but it hides differences between individual days.

# Calculate average hourly load by day
daily_hourly_load <- complete_data %>%
  mutate(
    hour = lubridate::hour(date_time)
  ) %>%
  group_by(date, hour) %>%
  summarise(
    avg_load = mean(load),
    .groups = "drop"
  ) %>%
  mutate(
    day_of_week = lubridate::wday(
      date,
      label = TRUE,
      abbr = FALSE
    )
  )

# Compare hourly load profiles across days
ggplot(
  daily_hourly_load,
  aes(
    x = hour,
    y = avg_load,
    group = date,
    color = day_of_week
  )
) +
  geom_line(linewidth = 1) +
  scale_x_continuous(
    breaks = 0:23,
    labels = sprintf("%02d:00", 0:23)
  ) +
  scale_color_manual(
    values = c(
      "Sunday" = "#ED6B5A",
      "Monday" = "#E88504",
      "Tuesday" = "#DDBB25",
      "Wednesday" = "#6EAD47",
      "Thursday" = "#49A3A9",
      "Friday" = "#AE85CE",
      "Saturday" = "#E884CD"
    )
  ) +
  labs(
    title = "Hourly BPA Load Profiles by Day",
    x = "Hour of Day",
    y = "Average Load (MW)",
    color = "Day"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(face = "bold")
  )

# Reshape data for time-series visualization
load_long <- load_data %>%
  pivot_longer(
    cols = c(load, ver, hydro, fossil_biomass, nuclear, interchange),
    names_to = "resource",
    values_to = "mw"
  ) %>%
  mutate(
    resource = recode(
      resource,
      load = "Load",
      ver = "VER",
      hydro = "Hydro",
      fossil_biomass = "Fossil/Biomass",
      nuclear = "Nuclear",
      interchange = "Net Interchange"
    )
  )

# Plot five-minute BPA load and generation data
ggplot(
  load_long,
  aes(
    x = date_time,
    y = mw,
    color = resource,
    group = resource
  )
) +
  geom_line(linewidth = 0.4, na.rm = TRUE) +

  scale_color_manual(
    values = c(
      "Load" = "#E41A1C",
      "VER" = "#228B22",
      "Hydro" = "#2878B5",
      "Fossil/Biomass" = "#8B6F47",
      "Nuclear" = "#7B3294",
      "Net Interchange" = "#777777"
    )
  ) +

  scale_x_datetime(
    date_breaks = "1 day",
    date_labels = "%b %d",
    timezone = "America/Los_Angeles",
    limits = range(load_data$date_time)
  ) +

  labs(
    title = "BPA Balancing Authority Load and Generation",
    subtitle = "Five-minute observations",
    x = "Date/Time",
    y = "MW",
    color = NULL
  ) +

  theme_minimal() +

  theme(
    plot.title = element_text(face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "top",
    panel.grid.minor = element_blank()
  )

# Define plot boundaries
plot_start <- min(load_data$date_time, na.rm = TRUE)
plot_end <- max(load_data$date_time, na.rm = TRUE)

# Set vertical position for weekend labels
label_y <- max(load_data$load, na.rm = TRUE) * 1.08

# Identify weekend days
weekend_days <- tibble(
  date = seq(
    as.Date(plot_start, tz = "America/Los_Angeles"),
    as.Date(plot_end, tz = "America/Los_Angeles"),
    by = "day"
  )
) %>%
  filter(lubridate::wday(date) %in% c(1, 7)) %>%
  mutate(
    # Define midnight boundaries in Pacific time
    xmin = as.POSIXct(
      paste(date, "00:00:00"),
      tz = "America/Los_Angeles"
    ),
    xmax = as.POSIXct(
      paste(date + 1, "00:00:00"),
      tz = "America/Los_Angeles"
    ),

    # Restrict shading to the visible plotting window
    visible_start = pmax(xmin, plot_start),
    visible_end = pmin(xmax, plot_end),

    # Center labels within the visible shaded region
    label_x = visible_start +
      (visible_end - visible_start) / 2,

    label = weekdays(date)
  )

# Check weekend dates and shading boundaries
weekend_days %>%
  select(date, xmin, xmax, visible_start, visible_end, label_x, label)

# Plot five-minute BPA load and generation data
ggplot() +

  # Weekend background shading
  geom_rect(
    data = weekend_days,
    aes(
      xmin = visible_start,
      xmax = visible_end,
      ymin = -Inf,
      ymax = Inf
    ),
    fill = "#D9D9D9",
    alpha = 0.45,
    inherit.aes = FALSE
  ) +

  # Five-minute load and generation observations
  geom_line(
    data = load_long,
    aes(
      x = date_time,
      y = mw,
      color = resource,
      group = resource
    ),
    linewidth = 0.4,
    na.rm = TRUE
  ) +

  # Weekend labels
  geom_text(
    data = weekend_days,
    aes(
      x = label_x,
      y = label_y,
      label = label
    ),
    inherit.aes = FALSE,
    fontface = "bold",
    size = 4
  ) +

  # Resource colors
  scale_color_manual(
    breaks = c(
      "Load",
      "VER",
      "Hydro",
      "Fossil/Biomass",
      "Nuclear",
      "Net Interchange"
    ),
    values = c(
      "Load" = "#E41A1C",
      "VER" = "#228B22",
      "Hydro" = "#2878B5",
      "Fossil/Biomass" = "#8B6F47",
      "Nuclear" = "#7B3294",
      "Net Interchange" = "#777777"
    )
  ) +

  # Date axis
  scale_x_datetime(
    date_breaks = "1 day",
    date_labels = "%a\n%b %d",
    timezone = "America/Los_Angeles",
    limits = c(plot_start, plot_end),
    expand = expansion(mult = 0)
  ) +

  # Y-axis
  scale_y_continuous(
    expand = expansion(mult = c(0.02, 0.12))
  ) +

  labs(
    title = "BPA Balancing Authority Load and Generation",
    subtitle = "Five-minute observations | November 7–14, 2026",
    x = "Date/Time (Pacific)",
    y = "Power (MW)",
    color = NULL,
    caption = "Source: BPA Balancing Authority Load & Total VER. Five-minute SCADA observations."
  ) +

  theme_minimal() +

  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 11),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "top",
    legend.box = "horizontal",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_line(
      color = "grey80",
      linewidth = 0.5
    ),
    plot.caption = element_text(
      color = "grey40",
      size = 8
    )
  )

Forecasting

In a typical regression we may estimate \[Load_t = \beta_0 + Temperature_t + \epsilon_t\] to investigate the relationship between load and temperature.

In forecasting, the objective would be to predict a future observation as a function of past observations: \[\widehat{Load_{t+1}} = f(Load_t, Load_{t-1}, Load_{t-24}, ...)\]

In order to practice, I will be using 2025 data as my training data for forecasting, and the 2026 data as testing data. Data is sourced from bullet 5 of https://transmission.bpa.gov/Business/Operations/Wind/.

# Import 2025 BPA data
data_25 <- read_excel(
  "~/Downloads/WindGenTotalLoadYTD_2025.xlsx",
  skip = 1
) %>%
  clean_names()

# Import 2026 BPA data
data_26 <- read_excel(
  "~/Downloads/WindGenTotalLoadYTD_2026.xlsx",
  skip = 1
) %>%
  clean_names()

# Compare missing values by column
missing_comparison <- tibble(
  variable = names(data_25),
  missing_2025 = colSums(is.na(data_25)),
  missing_2026 = colSums(is.na(data_26))
)

# Display missing-value comparison
missing_comparison %>%
  knitr::kable(
    caption = "Missing Values: 2025 vs. 2026",
    col.names = c("Variable", "2025 Missing", "2026 Missing")
  )

Missing Values: 2025 vs. 2026

Variable	2025 Missing	2026 Missing
date_time	0	0
date_time_utc	0	0
total_wind_base_schedule_forecast_in_bpa_control_area_mw_scada_187517	2	0
total_wind_basepoint_forecast_in_bpa_control_area_mw_scada_103349	2	0
total_wind_generation_in_bpa_control_area_mw_scada_79687	2	27
total_bpa_control_area_load_mw_scada_45583	2	0
total_hydro_generation_mw_scada_79682	2	27
total_fossil_biomass_generation_mw_scada_164377	2	27
total_nuclear_generation_mw_70681	2	0
net_interchange_mw_scada_45581	2	0
total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513	13438	9135
total_solar_basepoint_forecast_in_bpa_control_area_mw_scada_177165	2	0
total_solar_generation_in_bpa_control_area_mw_scada_142383	2	27

# Compare total observations
nrow(data_25)

#> [1] 105120

nrow(data_26)

#> [1] 91008

Combining and inspecting the raw data

# Combine 2025 and 2026 datasets
bpa_raw <- bind_rows(
  data_25 %>% mutate(source_year = 2025),
  data_26 %>% mutate(source_year = 2026)
) %>%
  arrange(date_time)


# Count duplicated timestamps
sum(duplicated(bpa_raw$date_time))

#> [1] 24

# Display any duplicated timestamps
bpa_raw %>%
  filter(duplicated(date_time) |
         duplicated(date_time, fromLast = TRUE))

# Parse local and UTC timestamps
bpa_raw <- bpa_raw %>%
  mutate(
    date_time = mdy_hm(
      date_time,
      tz = "America/Los_Angeles"
    ),
    date_time_utc = mdy_hm(
      date_time_utc,
      tz = "UTC"
    )
  )

bpa_raw %>%
  filter(
    as.Date(date_time_utc, tz = "UTC") == as.Date("2026-11-01")
  ) %>%
  select(date_time, date_time_utc) %>%
  print(n = 30, width = Inf)

#> # A tibble: 288 × 2
#>    date_time           date_time_utc      
#>    <dttm>              <dttm>             
#>  1 2026-10-31 17:00:00 2026-11-01 00:00:00
#>  2 2026-10-31 17:05:00 2026-11-01 00:05:00
#>  3 2026-10-31 17:10:00 2026-11-01 00:10:00
#>  4 2026-10-31 17:15:00 2026-11-01 00:15:00
#>  5 2026-10-31 17:20:00 2026-11-01 00:20:00
#>  6 2026-10-31 17:25:00 2026-11-01 00:25:00
#>  7 2026-10-31 17:30:00 2026-11-01 00:30:00
#>  8 2026-10-31 17:35:00 2026-11-01 00:35:00
#>  9 2026-10-31 17:40:00 2026-11-01 00:40:00
#> 10 2026-10-31 17:45:00 2026-11-01 00:45:00
#> 11 2026-10-31 17:50:00 2026-11-01 00:50:00
#> 12 2026-10-31 17:55:00 2026-11-01 00:55:00
#> 13 2026-10-31 18:00:00 2026-11-01 01:00:00
#> 14 2026-10-31 18:05:00 2026-11-01 01:05:00
#> 15 2026-10-31 18:10:00 2026-11-01 01:10:00
#> 16 2026-10-31 18:15:00 2026-11-01 01:15:00
#> 17 2026-10-31 18:20:00 2026-11-01 01:20:00
#> 18 2026-10-31 18:25:00 2026-11-01 01:25:00
#> 19 2026-10-31 18:30:00 2026-11-01 01:30:00
#> 20 2026-10-31 18:35:00 2026-11-01 01:35:00
#> 21 2026-10-31 18:40:00 2026-11-01 01:40:00
#> 22 2026-10-31 18:45:00 2026-11-01 01:45:00
#> 23 2026-10-31 18:50:00 2026-11-01 01:50:00
#> 24 2026-10-31 18:55:00 2026-11-01 01:55:00
#> 25 2026-10-31 19:00:00 2026-11-01 02:00:00
#> 26 2026-10-31 19:05:00 2026-11-01 02:05:00
#> 27 2026-10-31 19:10:00 2026-11-01 02:10:00
#> 28 2026-10-31 19:15:00 2026-11-01 02:15:00
#> 29 2026-10-31 19:20:00 2026-11-01 02:20:00
#> 30 2026-10-31 19:25:00 2026-11-01 02:25:00
#> # ℹ 258 more rows

# Calculate the difference between consecutive UTC timestamps
bpa_raw <- bpa_raw %>%
  arrange(date_time_utc) %>%
  mutate(
    interval_minutes = as.numeric(
      difftime(
        date_time_utc,
        lag(date_time_utc),
        units = "mins"
      )
    )
  )

# Count observations by interval length
bpa_raw %>%
  count(interval_minutes, sort = TRUE)

# Count missing values in each variable
missing_summary <- bpa_raw %>%
  summarise(
    across(
      everything(),
      ~ sum(is.na(.x))
    )
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "missing_values"
  ) %>%
  arrange(desc(missing_values))

knitr::kable(
  missing_summary,
  caption = "Missing Values by Variable"
)

Missing Values by Variable

variable	missing_values
total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513	22573
total_wind_generation_in_bpa_control_area_mw_scada_79687	29
total_hydro_generation_mw_scada_79682	29
total_fossil_biomass_generation_mw_scada_164377	29
total_solar_generation_in_bpa_control_area_mw_scada_142383	29
total_wind_base_schedule_forecast_in_bpa_control_area_mw_scada_187517	2
total_wind_basepoint_forecast_in_bpa_control_area_mw_scada_103349	2
total_bpa_control_area_load_mw_scada_45583	2
total_nuclear_generation_mw_70681	2
net_interchange_mw_scada_45581	2
total_solar_basepoint_forecast_in_bpa_control_area_mw_scada_177165	2
interval_minutes	1
date_time	0
date_time_utc	0
source_year	0

# Identify observations with missing measurements
missing_patterns <- bpa_raw %>%
  mutate(
    missing_count = rowSums(
      is.na(
        across(
          c(
            total_wind_generation_in_bpa_control_area_mw_scada_79687,
            total_bpa_control_area_load_mw_scada_45583,
            total_hydro_generation_mw_scada_79682,
            total_fossil_biomass_generation_mw_scada_164377,
            total_nuclear_generation_mw_70681,
            net_interchange_mw_scada_45581,
            total_solar_generation_in_bpa_control_area_mw_scada_142383,
            total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
          )
        )
      )
    )
  ) %>%
  filter(missing_count > 0) %>%
  select(date_time, source_year, missing_count)

# Summarize how many observations have each pattern
missing_patterns %>%
  count(source_year, missing_count)

# Missing solar base schedule forecast by year
bpa_raw %>%
  group_by(source_year) %>%
  summarise(
    total_observations = n(),
    missing_values = sum(is.na(
      total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
    )),
    percent_missing = mean(is.na(
      total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
    )) * 100,
    .groups = "drop"
  )

# Missing solar forecast by calendar date
bpa_raw %>%
  mutate(date = as.Date(date_time)) %>%
  group_by(source_year, date) %>%
  summarise(
    missing_values = sum(is.na(
      total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
    )),
    .groups = "drop"
  ) %>%
  filter(missing_values > 0) %>%
  arrange(desc(missing_values)) %>%
  head(20)

# Identify which variables are missing at each affected timestamp
bpa_raw %>%
  filter(
    if_any(
      c(
        total_wind_generation_in_bpa_control_area_mw_scada_79687,
        total_bpa_control_area_load_mw_scada_45583,
        total_hydro_generation_mw_scada_79682,
        total_fossil_biomass_generation_mw_scada_164377,
        total_nuclear_generation_mw_70681,
        net_interchange_mw_scada_45581,
        total_solar_generation_in_bpa_control_area_mw_scada_142383,
        total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
      ),
      is.na
    )
  ) %>%
  mutate(
    missing_variables = pmap_chr(
      across(
        c(
          total_wind_generation_in_bpa_control_area_mw_scada_79687,
          total_bpa_control_area_load_mw_scada_45583,
          total_hydro_generation_mw_scada_79682,
          total_fossil_biomass_generation_mw_scada_164377,
          total_nuclear_generation_mw_70681,
          net_interchange_mw_scada_45581,
          total_solar_generation_in_bpa_control_area_mw_scada_142383,
          total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
        )
      ),
      \(...) {
        values <- list(...)
        names(values)[map_lgl(values, is.na)] %>%
          paste(collapse = ", ")
      }
    )
  ) %>%
  count(source_year, missing_variables, sort = TRUE) %>%
  print(n = Inf, width = Inf)

#> # A tibble: 4 × 3
#>   source_year
#>         <dbl>
#> 1        2025
#> 2        2026
#> 3        2026
#> 4        2025
#>   missing_variables                                                             
#>   <chr>                                                                         
#> 1 total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513        
#> 2 total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513        
#> 3 total_wind_generation_in_bpa_control_area_mw_scada_79687, total_hydro_generat…
#> 4 total_wind_generation_in_bpa_control_area_mw_scada_79687, total_bpa_control_a…
#>       n
#>   <int>
#> 1 13436
#> 2  9135
#> 3    27
#> 4     2

# Examine solar forecast missingness by hour and month
solar_missing_hourly <- bpa_raw %>%
  mutate(
    year = year(date_time),
    month = month(date_time, label = TRUE, abbr = FALSE),
    hour = hour(date_time)
  ) %>%
  group_by(year, month, hour) %>%
  summarise(
    observations = n(),
    missing = sum(is.na(
      total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
    )),
    percent_missing = mean(is.na(
      total_solar_base_schedule_forecast_in_bpa_control_area_mw_scada_187513
    )) * 100,
    .groups = "drop"
  )

# Plot solar forecast missingness by hour
ggplot(
  solar_missing_hourly,
  aes(x = hour, y = month, fill = percent_missing)
) +
  geom_tile(color = "white", linewidth = 0.2) +
  facet_wrap(~ year) +
  scale_fill_gradient(
    low = "#FFF7EC",
    high = "#B30000",
    name = "% Missing"
  ) +
  scale_x_continuous(breaks = seq(0, 23, by = 2)) +
  labs(
    title = "Solar Base Schedule Forecast Missingness",
    subtitle = "Percentage of missing observations by hour and month",
    x = "Hour of Day",
    y = "Month"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(angle = 0)
  )

We need to convert the 5 minute into hourly load first.

# Create an hourly load dataset
hourly_load <- bpa_raw %>%
  mutate(
    hour_utc = floor_date(date_time_utc, unit = "hour")
  ) %>%
  group_by(hour_utc) %>%
  summarise(
    load_mw = mean(
      total_bpa_control_area_load_mw_scada_45583,
      na.rm = TRUE
    ),
    observations = sum(
      !is.na(total_bpa_control_area_load_mw_scada_45583)
    ),
    .groups = "drop"
  ) %>%
  # Keep hours with at least 11 of 12 five-minute load observations
  filter(observations >= 11) %>%
  mutate(
    date_time = with_tz(hour_utc, tzone = "America/Los_Angeles"),
    date = as.Date(date_time),
    year = year(date_time),
    month = factor(month(date_time, label = TRUE, abbr = FALSE)),
    day_of_week = factor(
      wday(date_time, label = TRUE, abbr = FALSE),
      levels = c(
        "Monday", "Tuesday", "Wednesday", "Thursday",
        "Friday", "Saturday", "Sunday"
      )
    ),
    hour = factor(hour(date_time))
  )

# Inspect the hourly dataset
glimpse(hourly_load)

#> Rows: 16,344
#> Columns: 9
#> $ hour_utc     <dttm> 2025-01-01 08:00:00, 2025-01-01 09:00:00, 2025-01-01 10:…
#> $ load_mw      <dbl> 6965.417, 6888.667, 6820.667, 6818.500, 6908.250, 7133.58…
#> $ observations <int> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 1…
#> $ date_time    <dttm> 2025-01-01 00:00:00, 2025-01-01 01:00:00, 2025-01-01 02:…
#> $ date         <date> 2025-01-01, 2025-01-01, 2025-01-01, 2025-01-01, 2025-01-…
#> $ year         <dbl> 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, 2025, 202…
#> $ month        <ord> January, January, January, January, January, January, Jan…
#> $ day_of_week  <ord> Wednesday, Wednesday, Wednesday, Wednesday, Wednesday, We…
#> $ hour         <fct> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…

# Check the number of hourly observations by year
hourly_load %>%
  count(year)

# Examine the first few observations
head(hourly_load, 10)

Model 1; Hourly BPA Load Forecasting

Our first model will use three kinds of information to predict hourly BPA load: 1. Yesterday’s load: What was demand at the same hour 24 hours ago? 2. Last week’s load: What was demand at the same hour seven days ago? 3. Calendar patterns: What hour, day of the week, and month is it? These are useful predictors because electricity demand has strong daily, weekly, and seasonal patterns.

# Create hourly load observations
hourly_data <- bpa_raw %>%
  mutate(
    hour_utc = floor_date(date_time_utc, unit = "hour")
  ) %>%
  group_by(hour_utc) %>%
  summarise(
    load_mw = mean(
      total_bpa_control_area_load_mw_scada_45583,
      na.rm = TRUE
    ),
    observations = sum(
      !is.na(total_bpa_control_area_load_mw_scada_45583)
    ),
    .groups = "drop"
  ) %>%
  mutate(
    date_time_local = with_tz(
      hour_utc,
      tzone = "America/Los_Angeles"
    ),
    hour = hour(date_time_local),
    weekday = wday(date_time_local, label = TRUE),
    month = month(date_time_local, label = TRUE)
  ) %>%
  arrange(hour_utc)

# Create lagged load variables
hourly_data <- hourly_data %>%
  mutate(
    load_lag_24 = lag(load_mw, 24),
    load_lag_168 = lag(load_mw, 168)
  )

# Define the testing period
test_start <- max(hourly_data$hour_utc) - days(30)

# Separate training and testing observations
training_data <- hourly_data %>%
  filter(
    hour_utc < test_start,
    !is.na(load_lag_24),
    !is.na(load_lag_168)
  )

testing_data <- hourly_data %>%
  filter(
    hour_utc >= test_start,
    !is.na(load_lag_24),
    !is.na(load_lag_168)
  )

# Check sample sizes
nrow(training_data)

#> [1] 15455

nrow(testing_data)

#> [1] 721

# Train hourly load forecasting model
load_model_1 <- lm(
  load_mw ~
    load_lag_24 +
    load_lag_168 +
    factor(hour) +
    factor(weekday) +
    factor(month),
  data = training_data
)

# Examine model results
summary(load_model_1)

#> 
#> Call:
#> lm(formula = load_mw ~ load_lag_24 + load_lag_168 + factor(hour) + 
#>     factor(weekday) + factor(month), data = training_data)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1093.12  -165.53    -4.74   156.59  1470.69 
#> 
#> Coefficients:
#>                     Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)        550.38941   32.33611  17.021  < 2e-16 ***
#> load_lag_24          0.81830    0.00469 174.490  < 2e-16 ***
#> load_lag_168         0.09443    0.00471  20.047  < 2e-16 ***
#> factor(hour)1      -11.23011   15.60992  -0.719 0.471892    
#> factor(hour)2      -15.80771   15.64396  -1.010 0.312288    
#> factor(hour)3      -16.01995   15.62963  -1.025 0.305392    
#> factor(hour)4       -4.44846   15.60536  -0.285 0.775602    
#> factor(hour)5       23.22122   15.65102   1.484 0.137912    
#> factor(hour)6       59.86187   15.93314   3.757 0.000173 ***
#> factor(hour)7       86.19769   16.29768   5.289 1.25e-07 ***
#> factor(hour)8       94.91714   16.46166   5.766 8.28e-09 ***
#> factor(hour)9       94.17757   16.46834   5.719 1.09e-08 ***
#> factor(hour)10      88.90377   16.38459   5.426 5.85e-08 ***
#> factor(hour)11      84.78254   16.31688   5.196 2.06e-07 ***
#> factor(hour)12      81.60688   16.26710   5.017 5.32e-07 ***
#> factor(hour)13      79.43408   16.23479   4.893 1.00e-06 ***
#> factor(hour)14      79.23895   16.22874   4.883 1.06e-06 ***
#> factor(hour)15      83.20548   16.29070   5.108 3.30e-07 ***
#> factor(hour)16      95.20178   16.49854   5.770 8.06e-09 ***
#> factor(hour)17     111.00478   16.79405   6.610 3.98e-11 ***
#> factor(hour)18     116.27646   16.87912   6.889 5.85e-12 ***
#> factor(hour)19     111.62842   16.79176   6.648 3.07e-11 ***
#> factor(hour)20     100.66294   16.58239   6.070 1.31e-09 ***
#> factor(hour)21      79.73851   16.22745   4.914 9.02e-07 ***
#> factor(hour)22      49.00788   15.84503   3.093 0.001985 ** 
#> factor(hour)23      19.71193   15.64367   1.260 0.207668    
#> factor(weekday).L -288.83999    6.15153 -46.954  < 2e-16 ***
#> factor(weekday).Q -174.26966    6.27317 -27.780  < 2e-16 ***
#> factor(weekday).C  115.91315    5.98498  19.367  < 2e-16 ***
#> factor(weekday)^4 -180.77676    6.02393 -30.010  < 2e-16 ***
#> factor(weekday)^5   78.43297    5.97850  13.119  < 2e-16 ***
#> factor(weekday)^6    4.07925    5.95945   0.685 0.493669    
#> factor(month).L    -12.94205    9.86977  -1.311 0.189782    
#> factor(month).Q    144.51632   10.75391  13.438  < 2e-16 ***
#> factor(month).C     -9.36014    8.24737  -1.135 0.256425    
#> factor(month)^4     65.06748    8.70285   7.477 8.04e-14 ***
#> factor(month)^5      4.62057    8.97494   0.515 0.606679    
#> factor(month)^6    -67.28172    8.67032  -7.760 9.02e-15 ***
#> factor(month)^7    -12.13790    8.15232  -1.489 0.136537    
#> factor(month)^8     42.02975    7.85768   5.349 8.98e-08 ***
#> factor(month)^9     21.65816    7.61716   2.843 0.004470 ** 
#> factor(month)^10    25.56769    7.38470   3.462 0.000537 ***
#> factor(month)^11     2.38185    7.32396   0.325 0.745024    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 280 on 15412 degrees of freedom
#> Multiple R-squared:  0.9128, Adjusted R-squared:  0.9126 
#> F-statistic:  3842 on 42 and 15412 DF,  p-value: < 2.2e-16

We are trying to forecast load_mw, or the hourly BPA balancing area electricity load, as measured in megawatts. We are using the information of load at the same hour of the previous day, same hour of the previous week, the hour of the day, the day of the week, and the month of the year. In other words, we are asking “given what electricity demand looked like recently, and knowing the time of day, day of the week, and time of year, how accurately can we predict electricity demand?

Before training the model we transformed the original 5-minute increments into hourly observations. This way we can work with an hour-interval lag as opposed to 5-minute lags.

We trained the model on all of the 2025 points, minus the final 30 days of observations. The test period is used for rolling one hour ahead predictions. At each hour the model can use actual load from the previous day and week. It is not forecasting all 30 days at once from a single starting point.

The model used is as follows: load_model_1 <- lm(load_mw ~ load_lag_24 + load_lag_168 + factor(hour) + factor(weekday) + factor(month), data = training_data)

Component	What it Captures
load_mw	The outcome we are predicting
load_lag_24	Previous day demand
load_lag_168	Previous week demand
factor(hour)	Differences between hours of the day
factor(weekday)	Differences between days of the week
factor(month)	Differences between months

Conceptually, this translates to \[Load_t = \beta_0 + \beta_1Load_{t-24} + \beta_2Load_{t-168} + Hour Effects + Weekday Effects + Month Effects + \epsilon_t\]

Interpreting Results

The previous-day load coefficient is 0.8183, and the previous-week load’s is 0.0944, both statistically significant. Interpreted, the previous-day coefficient means that, holding all other predictors constant, a 1 MW increase in yesterday’s load at the same hour is associated with an estimated 0.818 MW increase in predicted load. The model is relying heavily on the prior day’s load. The prior week contributes additional information, but the estimated coefficient is much smaller.

The hour-of-day effects use the hour 0 (midnight) as the reference category. The strongest hourly effect in this output is at 6pm, approximately 116 MW above the midnight reference. This is consistent with the idea that electricity demand has a daily pattern. However, these are adjustments relative to the reference hour, not the average load difference between those hours.

However, the model is treating the labelled weekday and monthly variables as ordered factors. Instead of estimating separate coefficient foe each weekday or month relative to a reference category, it is estimating polynomial contrasts. We can re-run this model using ordinary categorical contrasts by converting the factors into character first.

load_model_2 <- lm(
  load_mw ~
    load_lag_24 +
    load_lag_168 +
    factor(hour) +
    factor(as.character(weekday)) +
    factor(as.character(month)),
  data = training_data
)

summary(load_model_2)

#> 
#> Call:
#> lm(formula = load_mw ~ load_lag_24 + load_lag_168 + factor(hour) + 
#>     factor(as.character(weekday)) + factor(as.character(month)), 
#>     data = training_data)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1093.12  -165.53    -4.74   156.59  1470.69 
#> 
#> Coefficients:
#>                                    Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)                       393.47123   32.23982  12.205  < 2e-16 ***
#> load_lag_24                         0.81830    0.00469 174.490  < 2e-16 ***
#> load_lag_168                        0.09443    0.00471  20.047  < 2e-16 ***
#> factor(hour)1                     -11.23011   15.60992  -0.719 0.471892    
#> factor(hour)2                     -15.80771   15.64396  -1.010 0.312288    
#> factor(hour)3                     -16.01995   15.62963  -1.025 0.305392    
#> factor(hour)4                      -4.44846   15.60536  -0.285 0.775602    
#> factor(hour)5                      23.22122   15.65102   1.484 0.137912    
#> factor(hour)6                      59.86187   15.93314   3.757 0.000173 ***
#> factor(hour)7                      86.19769   16.29768   5.289 1.25e-07 ***
#> factor(hour)8                      94.91714   16.46166   5.766 8.28e-09 ***
#> factor(hour)9                      94.17757   16.46834   5.719 1.09e-08 ***
#> factor(hour)10                     88.90377   16.38459   5.426 5.85e-08 ***
#> factor(hour)11                     84.78254   16.31688   5.196 2.06e-07 ***
#> factor(hour)12                     81.60688   16.26710   5.017 5.32e-07 ***
#> factor(hour)13                     79.43408   16.23479   4.893 1.00e-06 ***
#> factor(hour)14                     79.23895   16.22874   4.883 1.06e-06 ***
#> factor(hour)15                     83.20548   16.29070   5.108 3.30e-07 ***
#> factor(hour)16                     95.20178   16.49854   5.770 8.06e-09 ***
#> factor(hour)17                    111.00478   16.79405   6.610 3.98e-11 ***
#> factor(hour)18                    116.27646   16.87912   6.889 5.85e-12 ***
#> factor(hour)19                    111.62842   16.79176   6.648 3.07e-11 ***
#> factor(hour)20                    100.66294   16.58239   6.070 1.31e-09 ***
#> factor(hour)21                     79.73851   16.22745   4.914 9.02e-07 ***
#> factor(hour)22                     49.00788   15.84503   3.093 0.001985 ** 
#> factor(hour)23                     19.71193   15.64367   1.260 0.207668    
#> factor(as.character(weekday))Mon  381.44711    8.65473  44.074  < 2e-16 ***
#> factor(as.character(weekday))Sat -156.96082    8.50648 -18.452  < 2e-16 ***
#> factor(as.character(weekday))Sun   58.79477    8.58580   6.848 7.78e-12 ***
#> factor(as.character(weekday))Thu   74.92730    8.44044   8.877  < 2e-16 ***
#> factor(as.character(weekday))Tue  193.16386    8.48212  22.773  < 2e-16 ***
#> factor(as.character(weekday))Wed   75.52602    8.45156   8.936  < 2e-16 ***
#> factor(as.character(month))Aug     45.13859   10.47619   4.309 1.65e-05 ***
#> factor(as.character(month))Dec    145.00957   13.27526  10.923  < 2e-16 ***
#> factor(as.character(month))Feb    105.31472   12.52439   8.409  < 2e-16 ***
#> factor(as.character(month))Jan    164.55224   12.86631  12.789  < 2e-16 ***
#> factor(as.character(month))Jul     77.02550   10.49074   7.342 2.21e-13 ***
#> factor(as.character(month))Jun     57.61167   10.45875   5.508 3.68e-08 ***
#> factor(as.character(month))Mar     57.08695   10.61246   5.379 7.59e-08 ***
#> factor(as.character(month))May     23.27254   10.40324   2.237 0.025297 *  
#> factor(as.character(month))Nov     95.74971   12.83087   7.462 8.95e-14 ***
#> factor(as.character(month))Oct     59.12605   11.45453   5.162 2.48e-07 ***
#> factor(as.character(month))Sep    -21.55200   10.50638  -2.051 0.040252 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 280 on 15412 degrees of freedom
#> Multiple R-squared:  0.9128, Adjusted R-squared:  0.9126 
#> F-statistic:  3842 on 42 and 15412 DF,  p-value: < 2.2e-16

With this fixed, we see adjusted differences after accounting for the lagged loads and other calendar variables. It is using Friday as a weekday reference, April as a monthly reference, and midnight (hour 0) as an hourly reference.

Again, the previous day load remains our strongest predictor. its coefficient is about 0.82, compared with the 0.09 for previous-week load. Both are highly statistically significant. Weekday differences are also substantial in this model. Monday’s estimated effect is approximately 381 MW above Friday, while Saturday is approximately 157 MW below Friday. This suggests meaningful differenced in demand patterns across the week.

We also see seasonal differences. January and December have positiove coefficients relative to April. This is consistent with seasonal variation in electricity demand, although we would need to investigate temperature and weather before attributing those differences to any particular cause.

The R-squared remained the same across models, alongside the residual standard error and the F statistic. Our model seems to fit the training data well, but we must test it to be sure. ### Testing the Model

We test for: - MAE (mean absolute error): measures the average absolute difference between actual and predicted values. - RMSE (root mean squared error): similar to MAE, but instead of taking the absolute value of each error, it squares the error, averages those squared errors, then takes the squared root. This essentially weights the larger errors more. - MAPE (mean absolute percentage error): expresses forecasting errors as percentages rather than MW.

# Generate predictions for the testing period
testing_data <- testing_data %>%
  mutate(
    predicted_load = predict(
      load_model_2,
      newdata = testing_data
    ),

    # Simple forecasting benchmarks
    naive_24 = load_lag_24,
    naive_168 = load_lag_168
  )

# Compare forecasting accuracy
model_accuracy <- tibble(
  model = c(
    "Regression model",
    "Previous-day baseline",
    "Previous-week baseline"
  ),

  MAE = c(
    mean(abs(testing_data$load_mw - testing_data$predicted_load)),
    mean(abs(testing_data$load_mw - testing_data$naive_24)),
    mean(abs(testing_data$load_mw - testing_data$naive_168))
  ),

  RMSE = c(
    sqrt(mean((testing_data$load_mw - testing_data$predicted_load)^2)),
    sqrt(mean((testing_data$load_mw - testing_data$naive_24)^2)),
    sqrt(mean((testing_data$load_mw - testing_data$naive_168)^2))
  ),

  MAPE = c(
    mean(abs((testing_data$load_mw - testing_data$predicted_load) /
      testing_data$load_mw)) * 100,

    mean(abs((testing_data$load_mw - testing_data$naive_24) /
      testing_data$load_mw)) * 100,

    mean(abs((testing_data$load_mw - testing_data$naive_168) /
      testing_data$load_mw)) * 100
  )
)

# Display results
knitr::kable(
  model_accuracy,
  digits = 2,
  caption = "Hourly Load Forecasting Accuracy"
)

Hourly Load Forecasting Accuracy

model	MAE	RMSE	MAPE
Regression model	192.40	238.42	2.74
Previous-day baseline	214.58	283.84	3.06
Previous-week baseline	308.89	390.02	4.44

Interpreting this, we see that our regression model including all of the variables were closer than simply relying on the previous day or week baseline. But, there is still some significant error in our predicting. One way to investigate is to calculate the forecast error for every observation and summarize it by hour of the day.

# Calculate forecasting errors
testing_data <- testing_data %>%
  mutate(
    forecast_error = load_mw - predicted_load,
    absolute_error = abs(forecast_error)
  )

# Summarize forecasting accuracy by hour
hourly_accuracy <- testing_data %>%
  group_by(hour) %>%
  summarise(
    MAE = mean(absolute_error),
    RMSE = sqrt(mean(forecast_error^2)),
    mean_error = mean(forecast_error),
    observations = n(),
    .groups = "drop"
  )

# Display results
knitr::kable(
  hourly_accuracy,
  digits = 2,
  caption = "Forecasting Accuracy by Hour of Day"
)

Forecasting Accuracy by Hour of Day

hour	MAE	RMSE	mean_error	observations
0	166.89	209.07	-6.36	30
1	177.19	218.46	-11.77	31
2	187.33	230.14	-10.49	30
3	197.99	244.35	-11.73	30
4	208.70	255.38	-18.71	30
5	234.76	277.47	-19.03	30
6	298.53	354.03	-10.89	30
7	332.34	383.67	-2.40	30
8	282.49	326.33	7.24	30
9	234.47	269.87	14.94	30
10	192.00	227.72	15.43	30
11	167.77	211.50	11.94	30
12	164.40	207.56	8.28	30
13	159.73	206.47	1.59	30
14	161.80	207.15	-5.35	30
15	170.64	213.58	-13.67	30
16	183.42	220.71	-11.17	30
17	178.56	216.09	-3.17	30
18	179.55	215.43	-2.57	30
19	183.44	213.09	-8.87	30
20	160.29	188.11	-7.47	30
21	140.80	167.33	-3.92	30
22	129.02	148.38	0.39	30
23	125.97	151.24	3.26	30

# Select the final seven days of the testing period
plot_start <- max(testing_data$hour_utc) - days(7)

forecast_plot_data <- testing_data %>%
  filter(hour_utc >= plot_start) %>%
  select(
    hour_utc,
    load_mw,
    predicted_load
  ) %>%
  pivot_longer(
    cols = c(load_mw, predicted_load),
    names_to = "series",
    values_to = "load"
  ) %>%
  mutate(
    series = recode(
      series,
      load_mw = "Actual Load",
      predicted_load = "Predicted Load"
    )
  )

# Plot actual and predicted load
ggplot(
  forecast_plot_data,
  aes(x = hour_utc, y = load, color = series)
) +
  geom_line(linewidth = 0.7) +
  scale_color_manual(
    values = c(
      "Actual Load" = "#0072B2",
      "Predicted Load" = "#D55E00"
    )
  ) +
  labs(
    title = "Actual vs. Predicted BPA Load",
    subtitle = "Final seven days of the testing period",
    x = "Date",
    y = "Load (MW)",
    color = NULL
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )