Forecasting HW1

Author

Teddy kelly

1 Select a Monthly Time Series

The time series I have chosen to analyze is the US Total Energy Monthly Production which I obtained from the U.S. Energy Information Administration website. Energy data exhibits large amounts of seasonality, so I think it would be interesting to analyze and generate forecasts for.

Also, since 2010, US energy production has increased much more compared to total energy consumption which has caused total energy exports to exceed total energy imports in recent years, resulting in the United States becoming a net exporter of total energy (“U.s. Energy Information Administration - EIA - Independent Statistics and Analysis,” n.d.). This is a very interesting trend, and for the judgemental forecasts section, I will analyze external factors that have led to the increase in total energy production and think about how the forecasts may change under new policies and changes to various other factors. Below, I have loaded the data, as well as all of the necessary packages for this project.

library(fpp3)
library(kableExtra)
library(patchwork)
rm(list=ls())
# Loading The Total Energy Monthly Production Time Series
energy <- read.csv("/Users/teddykelly/Downloads/MER_TE1.csv")

# Conveting the dataset into a tsibble and keeping only the total energy production data

energy_ts <- energy |>
  filter(Description == "Total Energy Production, Fossil Fuel Equivalent") |>
  mutate(year = YYYYMM %/% 100,
         month_num = YYYYMM %% 100) |>
  filter(month_num <= 12) |>
  mutate(month = yearmonth(sprintf("%d-%02d", year, month_num))) |>
  select(-year, -month_num) |>
  as_tsibble(index = month)

Important note:

Since the original data contained the annual total energy consumption for each year, I had to filter out those observations so that only the monthly observations were left. Also, the way the time variable was originally stored required me to convert it into a proper format that could be handled by the tsibble.
I utilized ChatGPT to do this and the YYYYMM %/% 100 command extracts the year component from the original date format and the YYYYMM %% 100 command extracts the month component.
The sprintf command combines the year and month components to obtain the proper year month format that yearmonth() can operate on.

2 Creating Training and Test Sets

There are 634 observations in the total energy production time series, so this corresponds to the first 507 observations belonging to the training set and the last 127 observations to the test set when performing an 80-20 split. The 507th observation is March 2015, so the observations from January 1973 to March 2015 will belong to the training set and the observations from April 2015 to October 2025 will belong to the testing set.

# Splitting the time series into 80% in the training set and 20% in the testing set
train <- energy_ts |> filter_index(~"Mar 2015")
test <- energy_ts |> filter_index("Apr 2015"~.)

Now, with the training and testing sets created, I will now explore and decompose the training set.

3 Explore and Decompose the Training Set

Below, I have plotted the original training set for US Total Energy production using autoplot().

train |> autoplot(Value) +
  labs(title = "Total Energy Production, Fossil Fuel Equivalent Training Set",
       subtitle = "Jan 1973 - Mar 2015",
       x = "Time (Months)",
       y = "Trillion Btus")

Analysis

Trend: There is a clear upward trend in the level of energy production from the beginning of the time series until the end, and specifically starting just before 2010.
Seasonality: There appears to be seasonality as there are regular patterns of cycles throughout the training set. There do appear to be some changes in the variance of the seasonality over time, suggesting that the time series may not be variance stationary and may require differencing.
Structural Breaks: There are a few structural breaks and especially towards the beginning of the training set. Those two big drops in energy production occurred during February of 1978 and April of 1981.

3.1 Decomposition

I will now decompose the training set that we see above into its trend-cycle, seasonality, and remainder components to understand how each component contributes to the time series more clearly. I will perform both classical and STL decomposition and compare the effectiveness of each method.

Classical Decomposition

For classical decomposition, I will use additive decomposition to simplify the analysis and interpretation. To do so, I have used the classical_decomposition() command within the model() function and specified the type to be additive.

classic_decomp <- train |> 
  model(classical_decomposition(Value, type = 'additive')) |>
  components() |> autoplot()

STL Decomposition

stl_decomp <- train |> 
  model(STL(Value ~ trend() + season())) |>
  components() |> autoplot()

Comparison

Comparison

Trend-cycle Components
- The trend-cycle components are largely the same, however, the STL trend component appears slightly smoother and the classical trend component is missing observations at the beginning and end of the time series because classical decomposition utilizes moving averages.
Seasonal Components
- The most obvious difference between the two decomposition methods is the seasonality component.
- For classical decomposition, the seasonal component remains constant and quite large over time. However, for STL decomposition, since I did not specify the window to be periodic, the seasonal component varies over time.
Remainder Components
- The remainder components are very similar, and largely exhibit random behavior throughout the training set apart from the beginning of the series with those two major structural breaks.
- However, the remainder component does appear to vary slightly more in the classical decomposition.

4 Apply Transformations to the Training Set

To determine whether or not a transformation needs to be applied to the training set, I have performed a unit root KPSS test to evaluate if the time series is stationary. If the training set is not stationary, then I will perform differencing on the total energy production values on the training set to prepare for the future forecasts. For a KPSS test, the null hypothesis is that the time series is stationary. So, if the p-value is greater than \(\alpha=0.05\), then we fail to reject the null hypothesis, and if the p-value is less than 0.05, then we can reject the null hypothesis and must perform a transformation so that the time series becomes stationary. Below is the code for performing the stationarity tests and differencing of the values.

# Kpss test on the training data
kpss_train <- train |> features(Value, unitroot_kpss)

# Adding the differenced value to the training set
train <- train |> mutate(diff_value = difference(Value))

# Testing stationarity of differenced data
kpss_train_diff <- train |> features(diff_value, unitroot_kpss)

# Applying differencing to the testing set
test <- test |> mutate(diff_value = difference(Value))

# Performing the kpss test on the testing set
kpss_test_diff <- test |> features(diff_value, unitroot_kpss)

# Combining the stationarity results into one table
kpss_table <- data.frame(
  "Series" = c("Original Train", "Differenced Train", "Differenced Test"),
  "kpss_stat" = c(kpss_train[[1]], kpss_train_diff[[1]], kpss_test_diff[[1]]),
  "p_value" = c(kpss_train[[2]], kpss_train_diff[[2]], kpss_test_diff[[2]])
)
kpss_table |> kable(digits = 2)

Series	kpss_stat	p_value
Original Train	5.10	0.01
Differenced Train	0.06	0.10
Differenced Test	0.05	0.10

The p-value for the original training set is 0.01 which is less than 0.05, indicating that we can reject the null hypothesis of stationary, meaning that the training set is not stationary and requires differencing.
The next row shows that differencing the values result in the training set being stationary since the p-value is 0.10 which is greater than 0.05, meaning that we fail to reject the null hypothesis of stationarity, and thus, the differenced training set is now stationary.
I also applied differencing to the total energy production values in the testing set, and the last row indicates that the differenced testing set is also stationary since the p-value is 0.01. Therefore, the methods I will fit on the training set can be used to generate forecasts for the testing set.

5 Model and Forecast Using the Training Set

5.1 Mean Method on Differenced Data (Drift Method)

For my baseline model, I want to utilize the Mean method on the differenced values for total energy production and then back-transform those forecasts to make projections on the original level of the energy production. Although I have already computed the differenced values, generating forecasts on the differenced data and then back-transforming to obtain the forecasts for the original level requires some additional steps which are not necessary.

Luckily, the method I just described is mathematically equivalent to applying the Drift method on the original values of the time series, and this can be implemented in just a few lines using the RW() function. Below are the first 5 observations of the fitted values obtained using the drift method against the original values in the training set.

# Fitting the drift method on the training set
fit_mean_diff <- train |>
  model(drift = RW(Value ~ drift()))

augment(fit_mean_diff) |> head(5) |> kable(digits = 2)

.model	month	Value	.fitted	.resid	.innov
drift	1973 Jan	5404.72	NA	NA	NA
drift	1973 Feb	5155.12	5409.04	-253.92	-253.92
drift	1973 Mar	5419.56	5159.44	260.12	260.12
drift	1973 Apr	5160.81	5423.88	-263.07	-263.07
drift	1973 May	5411.25	5165.14	246.11	246.11

We can see that there is no fitted value for the first observation because fitted values for the training set are essentially one step ahead predictions using the previous observation. Since for the first observation, there is no previous value, then the drift method cannot generate a fitted value for the first observation.
Just looking at the pattern of the first 5 fitted values, it appears they are oscillating to match the level of the previous observation, indicating that the fitted values should closely align with the seasonality seen in the training set.
In fact, to visually observe this pattern, I have graphed the fitted values against the actual total energy production data in the training set below:

# Mean Difference Method
augment(fit_mean_diff) |>
    ggplot(aes(x = month)) +
    geom_line(aes(y = Value, color = "Actual Data")) +
    geom_line(aes(y = .fitted, color = "Fitted Values"), linetype = 'dashed') +
    scale_color_manual(
        values = c("black", "blue"),
        breaks = c("Actual Data", "Fitted Values")
    ) +
  theme(legend.position = 'bottom') +
  labs(title = "Drift Method Fitted Values Vs Actual Training Set Values",
       subtitle = "Us Total Energy Production from Jan 1973 - Mar 2015",
       x = "Time (Months)",
       y = "Total Energy Production in BTUs")

Looking at the graph above, we can see that even for the simple drift method, which is just the mean method applied to the differenced values and then back-transformed to the original level, it very accurately matches the trend and seasonality of actual total energy production in the training set.
The drift method does not actually account for seasonality, however it only appears to “capture” the seasonality in the training set because it uses the past observed value to make a one step ahead prediction.
However, when generating forecasts for unobserved future values, the drift method will no longer be able to “capture” the seasonality, but only the trend since it will no longer be able to use past observed values to make predictions..

5.2 Regression Model with Trend and Seasonal Dummy

The next model I will fit on the training data to generate forecasts for the testing set is a linear regression model that contains trend and seasonal dummy variables. Since the time series I am analyzing only contains the forecast variable that I am interested in studying (US Total Energy Production), the predictor variables in this case will be the trend dummy and the seasonal dummies. Since this is a monthly time series, R will automatically generate 11 seasonal dummy variables with one of the months being captured by the intercept term. Whichever month is captured by the intercept term will be the baseline month that we will compare the coefficient values for the other monthly dummy variables to.

# Fitting the regression model with trend and season dummies on training set
fit_regression <- train |>
  model(TSLM(diff_value ~ trend() + season()))

# Backtransforming using the lag of the orginial value

# Displaying the regression table output
fit_regression |> report()

Series: diff_value 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-805.873  -94.630   -5.589   88.953  612.529 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      79.14734   26.12256   3.030  0.00258 ** 
trend()           0.02262    0.04589   0.493  0.62226    
season()year2  -555.63141   32.70664 -16.988  < 2e-16 ***
season()year3   412.05583   32.70635  12.599  < 2e-16 ***
season()year4  -391.99138   32.90025 -11.915  < 2e-16 ***
season()year5    89.54759   32.89971   2.722  0.00672 ** 
season()year6  -183.62472   32.89923  -5.581 3.94e-08 ***
season()year7   -35.83396   32.89881  -1.089  0.27659    
season()year8    70.39480   32.89846   2.140  0.03286 *  
season()year9  -350.70475   32.89817 -10.660  < 2e-16 ***
season()year10  103.30958   32.89795   3.140  0.00179 ** 
season()year11 -226.91549   32.89779  -6.898 1.63e-11 ***
season()year12  102.00907   32.89769   3.101  0.00204 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 150.8 on 493 degrees of freedom
Multiple R-squared: 0.752,  Adjusted R-squared: 0.746
F-statistic: 124.6 on 12 and 493 DF, p-value: < 2.22e-16

The adjusted R-squared value is about 0.75 meaning that about 75% of the variation in US total energy consumption is explained by the time series linear regression model which is fairly high.
The F-statistic is also high with a value of about 125 and is statistically significant, meaning that there is evidence of a relationship between the response variable, US total energy production, and its predictor variables, which in this case are the trend and monthly dummy variables

Coefficient Interpretations

The table above displays coefficient estimates for each of the trend and monthly dummy variables. For the monthly dummies, R left out season()year1 which corresponds with January, so January will be the baseline month that is captured by the intercept term. For each monthly dummy variable, the coefficient represents how much US total energy production changed from the previous month on average compared to January since the model was fit on the differenced values. Below, I have provided coefficient interpretations for the intercept, trend, and a positive and negative monthly dummy variable:

Intercept (January)
- The intercept term corresponds to January, so the coefficient represents the predicted average monthly change in US energy production during January.
- The data measures total energy production in trillions of BTUs, so the coefficient \(\beta_0=79.15\) means that on average, from December to January, US total energy production increases by about 79 trillion BTUs.
- This coefficient is statistically significant since its p-value is below \(\alpha=0.05\), meaning that it’s unlikely that the coefficient estimate was obtained if the true population parameter was zero.
Trend
- The trend coefficient is positive indicating an upward trend in total energy production, and specifically, the coefficient of \(\beta_1=0.02\) means that on average, US total energy production increases by about 0.02 trillion BTUs every year.
- This coefficient is not statistically significant, however, when looking at the graph of the training set, it does appear that there is an upward trend.
- Also, I did experiment with specifying the linear regression without the trend dummy, and the forecasts obtained with that model did not capture the trend of the testing set as well as the regression model with the trend dummy.
season()year2 (February):
- The coefficient estimate is negative, meaning that February typically sees a lower change in total energy production compared to January. Specifically, the coefficient estimate of \(\beta_2=-555.63\) means that on average, total energy production changes by about 555 trillion BTUs less from January to February than it does from December to January.
- This coefficient is statistically significant, meaning that there is evidence of a difference in the change of total energy production during February compared to January.
season()year3 (March):
- The coefficient estimate is positive, meaning that March typically sees a higher change in total energy production compared to January. Specifically, the estimate of \(\beta_3=412\) means that on average, total energy production changes by about 412 trillion BTUs more from February to March than it does from December to January.
- The coefficient is also statistically significant suggesting that there is evidence of a difference in the change of total energy production during March compared to January.

Below, I have graphed the fitted values from the time series regression against the actual values from the training set.

actual_train <- train$Value

reg_aug <- augment(fit_regression) |> 
  mutate(actual_values = actual_train) |> 
  select(month, .fitted, actual_values)

reg_aug <- reg_aug |> mutate(yesterday_actual = lag(actual_values), original_fitted = yesterday_actual + .fitted)


reg_aug |>
    ggplot(aes(x = month)) +
    geom_line(aes(y = actual_values, color = "Actual Data")) +
    geom_line(aes(y = original_fitted, color = "Fitted Values"), linetype = 'dashed') +
    scale_color_manual(
        values = c("black", "blue"),
        breaks = c("Actual Data", "Fitted Values")
    ) +
  theme(legend.position = 'bottom') +
  labs(title = "Linear Regression Fitted Values Vs Actual Training Set Values",
       subtitle = "Us Total Energy Production from Jan 1973 - Mar 2015",
       x = "Time (Months)",
       y = "Total Energy Production in BTUs")

Looking at the graph above, it appears that the time series regression fitted values capture the trend and seasonality very well in the actual values from the training set. The fitted values from the regression model are very similar to the fitted values from the drift model, indicating that they both capture the trend and seasonality well in the training set.
However this does not mean that they will produce similar forecasts on the testing set, and in the next section, we will compare their forecasting accuracy and select the best model of the two

6 Evaluate Performance on Test Set

Below, I have generated the forecasts for both the drift and regression model and graphed their forecasts. This first graph includes both the training and testing set, and the second graph compares the forecasts only against the testing set.

# Applying forecasts the length of the test set for the drift method
mean_diff_fc <- fit_mean_diff |>
  forecast(h = nrow(test))


# Applying forecasts for the time series regression
regression_fc <- fit_regression |>
  forecast(h = nrow(test))

# Back transforming the forecasts on the differenced values to the original values
last_value <- last(train$Value)
regression_fc <- regression_fc |>
  mutate(original_forecasts = last_value + cumsum(.mean))

# Graphing the forecasrts from both methods against the train and test sets
train |>
  ggplot(aes(x = month)) +
  geom_line(aes(y = Value, color = "Training Set")) +
  geom_line(data = regression_fc, aes(y = original_forecasts, color = "Regression Forecasts")) +
  geom_line(data = mean_diff_fc, aes(y = .mean, color = "Drift Forecasts")) +
  geom_line(data = test, aes(y = Value, color = "Testing Set")) +
  scale_color_manual(
    values = c('black', 'blue', 'orange', 'red'),
    breaks = c("Training Set", "Regression Forecasts", "Drift Forecasts", "Testing Set")
  ) +
  theme(legend.position = 'bottom') +
  labs(title = "Drift and Time Series Regression Forecasts Against Training and Testing Sets",
       subtitle = "US Total Energy Production from Jan 1973 - Oct 2025",
       x = "time (Months)",
       y = "Total Energy Production (Trillions of BTUs")

# Graphing the forecasts against the testing set
#mean_diff_fc |> autoplot(train) +
 # autolayer(test, Value, color = 'red')

# Graphing the forecasts on just the testing set
test |>
  ggplot(aes(x = month)) +
  geom_line(aes(y = Value, color = "Testing Set")) +
  geom_line(data = regression_fc, aes(y = original_forecasts, color = "Regression Forecasts")) +
  geom_line(data = mean_diff_fc, aes(y = .mean, color = "Drift Forecasts")) +
  scale_color_manual(
    values = c('black', 'blue', 'red'),
    breaks = c("Testing Set", "Regression Forecasts", "Drift Forecasts")
  ) +
  theme(legend.position = 'bottom') +
  labs(title = "Drift and Time Series Regression Forecasts Against Testing Set Only",
       subtitle = "US Total Energy Production from Mar 2015 - Oct 2025",
       x = "time (Months)",
       y = "Total Energy Production (Trillions of BTUs")

Observing the graphs above, the forecasts produced from the linear regression model appear to be more accurate than those produced by the drift method. I was unable to include the prediction intervals in these particular plots, however when I graph the forecasts individually, the prediction intervals are much narrower for the linear regression model than the drift model, further suggesting the enhanced accuracy from the regression model.
Not only is the regression model able to match the seasonality of the actual data, it also better captures the trend compared to the drift method. In the next section, we will compare the accuracy metrics between the two methods to verify what we see in the graphs.
As mentioned in section 4, the drift method does not take into account seasonality, so for the forecasts, it only is able to capture some of the trend.
Although the regression model is not able to fully capture the trend of the actual data since it assumes a linear trend, we can see that the regression model accurately captures the seasonality of the actual values. Looking at the second graph in particular, whenever the actual values in the testing set increase, the forecasted values produced from the regression model also increase, and whenever the actual values decrease, the regression model correctly forecasts a decrease as well.

6.1 Comparing Accuracy Metrics

Below is a table comparing the accuracy metrics between the forecasts obtained from the drift method and the linear regression model. Since I had to manually back-transformed the forecasted differenced values from the regression model back to the original level of the series, I also had to manually compute the accuracy metrics for the regression forecasts.

mean_metrics <- accuracy(mean_diff_fc, energy_ts) |> select(.model, RMSE, MAE, MAPE)

# Comparing back-transformed forecasts to original data
original_data <- regression_fc |> as_tibble() |> select(month, original_forecasts) |>
    left_join(energy_ts |> as_tibble() |> select(month, Value), by = "month")

# Manually computing the accuracy metrics
regression_metrics <- original_data |> mutate(
    error = Value - original_forecasts,
    abs_error = abs(error),
    pct_error = abs(error / Value)
) |>
    summarize(
      .model = "Regression Model",
        RMSE = sqrt(mean(error**2, na.rm = T)),
        MAE = mean(abs_error, na.rm = T),
        MAPE = 100 * mean(pct_error, na.rm = T)
    )

# Bind the metrics together
metrics <- rbind(mean_metrics, regression_metrics)
colnames(metrics)[1] <- "Models"
metrics[[1,1]] <- "Drift Method"

metrics |> kable(digits = 2)

Models	RMSE	MAE	MAPE
Drift Method	730.99	622.54	7.37
Regression Model	422.98	365.18	4.44

As seen in the table above, the model accuracy metrics confirm what we found in the figures above that the time series regression model is able to forecast the actual values of the total energy production much more effectively.
In fact, the regression model has a significantly lower root mean square error and mean absolute error, and a lower mean absolute percentage error.
Therefore, between these two models, the time series linear regression is better for forecasting future US total energy production since it produces forecasts with less error than the drift method.

7 Judgmental Forecasts

Up to this point, we have analyzed the US total energy production time series and compared the forecasting accuracy of the drift method and a time series linear regression model with trend and seasonal dummy variables. Now, it’s now time to think about how these forecasts could potentially change if the circumstances were different. Different circumstances could mean a new policy that was implemented which could affect energy production or an economic shock that reduces the supply of the factors of production that are used to produce energy.

To do so, we will use scenario-based forecasting.

7.1 Scenario-Based Forecasts

There are three scenarios that I will consider for determining how the forecasts would change: A baseline case, an optimistic case, and a pessimistic case.

Optimistic Case:

For the optimistic case, US total energy production would continue to increase. Increases in total energy consumption is desired because it signals economic growth and means that there will be a greater supply of energy for American citizens to use. Also, I mentioned previously how since total energy production began significantly increasing in 2010, the US has now become a net exporter of total energy which is another positive effect since it will result in a trade surplus and increase GDP since net exports are rising.
This optimistic case of total energy production increasing could occur under a couple of circumstances. First, there could be technological innovation that makes extracting energy easier and more efficient, or there could be a rise in availability of renewable energy sources. This would be particularly beneficial for the US because not only is total energy production increasing which will lead to economic growth, but also the renewable energy would preserve the environment and allow long-term growth for future generations.
Another cause could be from an increase in total energy demand for consumption because the more people need to consume energy, the more need there will be for an increase in energy production. If we knew that the population was going to continue to increase and that more and more people will be living lifestyles that require large energy usage, then we would adjust our models to account for this information so that they could generate forecasts of increasing total energy production.
Specifically, for either technological innovation or an increase in energy demand, the drift method would need to be adjusted so that its slope would be steeper and capture this increase in production. Also, the regression model would likely have to be modified to be maybe piecewise linear so that it could change slopes to account for the increase in energy production.

Baseline Case

In the baseline case, the recent increase in US total energy production could level off and future total energy consumption would stay about the same. This could result from a reduction in total energy demand that causes the rate of total energy production increase to slow down since there isn’t a need to produce as much.
Also, perhaps there are no new technological innovations made that increases the availability of energy or lowers the cost of production.
To account for these circumstances, we would have to modify our forecasting methods so that their trends are close to non-existent. The drift method in particular would likely not work since if there is any type of trend in the training set, it will continue that trend for future unobserved values. Hence, applying the mean method on the original values could yield some fairly accurate forecasts if there is not much of a trend in the future.
For the regression method, again adjusting it to be piecewise linear would allow it to change slopes when the total energy production levels off. Also, removing the trend component from the model could reduce the positive trend that picked up from the training set and provide more accurate forecasts.

Pessimistic Case

In the pessimistic case, the US total energy production would decrease drastically. This would be devastating because it would indicate a potential economic downturn and it would mean that there would be a decrease in the supply of energy for US citizens. Also, this would mean that exports of total energy would decrease which could result in a trade deficit for total energy.
This decrease in total energy production could result from policy changes regarding energy in the United States. For example, President Trump has targeted renewable energy initiatives which could have a major impact on renewable energy production. Now, the US does not rely heavily on renewable energy, so total energy production will not likely decrease because of this. However, the decrease in renewable energy production could still have negative effects especially for the future.
Another cause of a decrease in US total energy could be an economic recession. If we look at the graph with the forecasts against the testing set, we can see that there was a significant drop in total energy production around the time of the pandemic. Energy production was able to recover after this, however if there is evidence suggesting a long period of economic turmoil, we might adjust our models to forecasts for a drop in energy production.
Adjusting the regression model to be non-linear so it can adjust its slope over time could be a potential solution.

In summary, judgmental forecasts are important in forecasting because there are many different scenarios that can occur in the future and there is no way of knowing which scenario will occur with complete certainty. Scenario-based forecasting allows us to account for these different possibilities and analyze how forecasted values would change depending on which scenario becomes reality. Based on these judgmental forecasts, policy makers can develop strategies accordingly and be prepared for any circumstance.

8 Interpretive Summary

In this study, I analyzed the US total energy production time series from January of 1973 to October of 2025. I implemented two forecasting methods (drift and linear regression model) which I fitted on the training set and then used to forecast future energy production values against the testing set.
Since the original time series was not stationary, I differenced the data to achieve stationary before fitting my forecasting models. Since I originally wanted to use the Mean method, I ended up using the drift method on the original values because this is mathematically equivalent to applying the mean method on the differenced data and then back-transforming the fitted and forecast value to the original level of the time series.
For the linear regression model, I used a trend dummy and monthly dummies as the predictors and I fit the model on the differenced values of the data. Overall, the linear regression model had an adjusted R-squared value of about 0.75 and a statistically significant F-statistic value of 125, indicating that the model fit the data well. I then manually back-transformed the fitted and forecasted values to their original levels and compared those forecasts to the forecasts produced from the drift method.
Overall, the regression model produced much more accurate forecasts of energy production than the drift method. The regression model had lower values for the RMSE, MAE, and MAPE.
Implementing judgement is very important when it comes to forecasting because it requires the forecaster to think about what else is going on in the world and how circumstances could affect forecasts. For the drift and regression models, the forecasts that were produced were fairly accurate in that they both captured the upward trend in total energy production in the US. However, going further into the future, as the models are specified now, they would not be able to account for any scenarios in which circumstances cause changes to US total energy production growth.

References

“U.s. Energy Information Administration - EIA - Independent Statistics and Analysis.” n.d. https://www.eia.gov/totalenergy/data/browser/index.php?tbl=TE1#/?f=M&start=197301&end=202510&charted=4-6-7-14.