Time Series Analysis of New York Harbor Conventional Gasoline Regular Spot Price

STATS 531 FINAL PROJECT

1. Introduction

Gasoline prices play a critical role in both global and local economies, influencing consumer spending, transportation costs, and broader macroeconomic trends. Because gasoline is a key input for many industries, its price fluctuations can serve as an economic indicator, with sharp increases often contributing to inflation. Given its economic significance, governments tightly regulate gasoline prices to mitigate volatility. For example, during recent economic downturns, the U.S. released one million barrels from its gasoline reserve to curb price surges ^[1], while OPEC+ cut oil production to raise prices. Due to such heavy government intervention^[2], we hypothesize that the leverage effect — the inverse relationship between asset returns and their volatility — is more limited in gasoline prices compared to freely traded commodities or financial assets.

This report analyzes historical monthly spot prices for conventional regular gasoline at New York Harbor^[3][4], covering the period from June 1986 to March 2025. To evaluate our hypothesis, we estimate and compare three distinct models:

• A modified stochastic volatility (SV) model with leverage effects

• A modified basic stochastic volatility model without leverage effects

• A GARCH model, which serves as the benchmark for comparison.

By assessing these models, we aim to determine whether the leverage effect—typically observed in freely traded assets—is diminished in gasoline prices due to regulatory influences.

2. Results

2.1 Data Preprocessing

The spot prices for conventional regular gasoline at New York Harbor remained relatively stable from 1986 to 1999 before entering a period of exponential growth. This upward trend was interrupted by the 2008 financial crisis, after which prices recovered and surged again. Between 2012 and 2013, multiple factors—including strong economic growth in the U.S. and China, as well as accommodative monetary policies in Europe—drove global crude oil demand upward, sustaining high petroleum prices. However, beginning in mid-2014, prices plummeted due to a global supply glut, with the decline accelerating through 2015^[5]. In 2020, gasoline prices experienced another sharp drop as the COVID-19 pandemic crushed demand. Yet the market rebounded quickly, and by early 2023, prices had risen sharply, contributing significantly to the inflationary pressures of 2022-2023 (Figure 1).

Let \(\{z_n, n = 1,...,N\}\) denote the gasoline price index series. We compute the demeaned log returns as:

\[ y^*_n = \log(z_n) - \log(z_{n-1}) \]

Figure 2 presents the log returns for both daily and monthly gasoline prices at New York Harbor from June 1986 to March 2025. While the daily returns appear smoother due to higher sampling frequency, our analysis focuses on monthly prices, where black swan events during the 2008 recession (time index 275) and 2020 pandemic (time index 405) periods become more pronounced. These features will prove particularly relevant when we adapt Breto's stochastic volatility model to account for such regime changes in gasoline price dynamics.

Note: All POMP framework instances described below were executed at run level 3.

2.2 Breto's (2014) SV - Leverage Model (POMP Framework)

2.2.1 Formulas

Here, volatility is modeled as a latent stochastic process, partially observed via the returns. We present a pomp implementation of Breto (2014)^[6][7] as

\[\begin{align} R_n &= \frac{\exp\{G_n\}-1}{\exp\{G_n\}+1} & (1)\\ G_n &= G_{n-1}+\nu_n & (2)\\ H_n &= \mu_h(1-\phi)+\phi H_{n-1}+\beta_{n-1}R_n\exp(-H_{n-1}/2)+\omega_n & (3)\\ Y_n &= \exp\left\{\frac{H_n}{2}\right\}\sigma_n & (4) \end{align}\]

where \(Y_n\) is the observed return, \(\beta_n=Y_n\sigma_\eta \sqrt{1-\phi^2}\), \(\{\epsilon_n\}\) is an i.i.d. \(N(0,1)\) sequence, \(\{\nu_n\}\) is an i.i.d. \(N(0,\sigma_{\nu}^2)\) sequence and \(\{\omega_n\}\) is \(N(0,\sigma_{\omega,n}^2)\) sequence where \(\sigma_{\omega,n}^2=\sigma_\eta^2(1-\phi^2)(1-R_n^2)\). \(H_n\) is the log volatility, \(G_n\) is Gaussian random walk. \(R_n\) is the leverage, defined on day \(n\) as the correlation between index return on day \(n-1\) and the increase in the log volatility from day \(n-1\) to day \(n\).

2.2.2 Simulation

# Set up pomp model using Breton et al. (2014) model
N_breto_statenames <- c("H","G","Y_state")
N_breto_rp_names <- c("sigma_nu","mu_h","phi","sigma_eta")
N_breto_ivp_names <- c("G_0","H_0")
N_breto_paramnames <- c(N_breto_rp_names, N_breto_ivp_names)

N_breto_rproc1 <- "
  double beta,omega,nu_noise;
  omega = rnorm(0,sigma_eta * sqrt( 1- phi*phi ) * 
    sqrt(1-tanh(G)*tanh(G)));
  nu_noise = rnorm(0, sigma_nu); // Change this
  G += nu_noise; // Change This
  beta = Y_state * sigma_eta * sqrt( 1- phi*phi );
  H = mu_h*(1 - phi) + phi*H + beta * tanh( G ) 
    * exp(-H/2) + omega;
"
N_breto_rproc2.sim <- "
  Y_state = rnorm( 0,exp(H/2) );
 "

N_breto_rproc2.filt <- "
  Y_state = covaryt;
 "

N_breto_rproc.sim <- paste(N_breto_rproc1, N_breto_rproc2.sim)
N_breto_rproc.filt <- paste(N_breto_rproc1, N_breto_rproc2.filt)

N_breto_rinit <- "
  G = G_0;
  H = H_0;
  Y_state = rnorm( 0,exp(H/2) );
"

N_breto_rmeasure <- "
   y=Y_state;
"

N_breto_dmeasure <- "
   lik=dnorm(y,0,exp(H/2),give_log);
"

N_breto_partrans <- parameter_trans(
  log=c("sigma_eta","sigma_nu"),
  logit="phi"
)

N_breto_filt <- pomp(data=data.frame(
    y=new_york_gasoline.ret.demeaned,time=1:length(new_york_gasoline.ret.demeaned)),
  statenames=N_breto_statenames,
  paramnames=N_breto_paramnames,
  times="time",
  t0=0,
  covar=covariate_table(
    time=0:length(new_york_gasoline.ret.demeaned),
    covaryt=c(0,new_york_gasoline.ret.demeaned),
    times="time"),
  rmeasure=Csnippet(N_breto_rmeasure),
  dmeasure=Csnippet(N_breto_dmeasure),
  rprocess=discrete_time(step.fun=Csnippet(N_breto_rproc.filt),
    delta.t=1),
  rinit=Csnippet(N_breto_rinit),
  partrans=N_breto_partrans
)

N_breto_params_test <- c(
  sigma_nu = exp(-4.5),  
  mu_h = -5.5,
  phi = expit(-0.5),     
  sigma_eta = exp(-10),
  G_0 = 0,
  H_0 = 0
)
  
N_breto_sim1.sim <- pomp(N_breto_filt, 
  statenames=N_breto_statenames,
  paramnames=N_breto_paramnames,
  rprocess=discrete_time(
    step.fun=Csnippet(N_breto_rproc.sim),delta.t=1)
)

N_breto_sim1.sim <- simulate(N_breto_sim1.sim,seed=1,params=N_breto_params_test)

N_breto_sim1.filt <- pomp(N_breto_sim1.sim, 
  covar=covariate_table(
    time=c(timezero(N_breto_sim1.sim),time(N_breto_sim1.sim)),
    covaryt=c(obs(N_breto_sim1.sim),NA),
    times="time"),
  statenames=N_breto_statenames,
  paramnames=N_breto_paramnames,
  rprocess=discrete_time(
    step.fun=Csnippet(N_breto_rproc.filt),delta.t=1)
)

We conducted parameter simulation studies to initialize our model fitting procedure. The initial parameter configuration \[ \theta_0 = (\sigma_\nu,\mu_h,\phi,\sigma_\eta,G_0,H_0) = (\exp(4.5), -5.5, \textrm{expit}(-0.5), \exp(-10),0,0) \] produced simulated trajectories that closely matched the observed data, except at structural break points corresponding to the 2008 recession (t=275) and 2020 pandemic (t=405) (Figure 3).

Using this initialization, we computed the filtered log-likelihood as 410.657 (SE = 0.079) for the simulated data. We subsequently performed local optimization around these parameter values to refine our estimates.

2.2.3 Local Search

Figure 4 and 5 show that our local search achieves reasonable convergence in log-likelihood, with differences between searches typically within a few log units. While parameters \(\sigma_\eta, \mu_h\), and \(\phi\) converge well, the pair plot suggests no linear relationship exists between any parameter and the log-likelihood. Notably, all searches yield an effective sample size (ESS) below 10 at the final iteration (time index 405, corresponding to the COVID-19 pandemic), indicating potential model misspecification.

This breakdown likely stems from the model's inability to capture extreme “black swan” events like the pandemic-induced price shock — a limitation also observed in our simulations. The issue is exacerbated by our use of monthly (rather than daily) returns, which amplify black swan events due to lower data granularity. While daily returns would provide smoother demeaned data and reduce outlier effects, our monthly sampling makes sudden price shifts more pronounced, complicating model fitting.

Given these findings, model adjustments are warranted. However, before modifying the specification, we first test whether a more extensive global search—with additional search rounds—might improve performance under these extreme conditions.

2.2.4 Global Search

Despite employing a global search (Figure 6 and 7), we observe the same limitations present in our local search results. The estimation procedure achieved consistent convergence, with the log-likelihood progressively improving across iterations to reach a maximum value of 429.8. However, we observe concerningly low effective sample sizes (ESS < 10) during the pandemic period in the final iteration, indicating potential difficulties in properly characterizing this extreme event. This empirical evidence motivates necessary model adjustments to better account for tail-risk scenarios like the COVID-19 market shock.

2.3 Modified SV - Leverage Model (POMP Framework)

2.3.1 Formulas

We retain the core structure from equations (1)-(3) but propose two key modifications to better capture extreme events and regime shifts:

Heavy-Tailed Observation Process^[8]

The observed return is modeled as: \[\begin{align} Y_n &= \exp\left\{\frac{H_n}{2}\right\}\sigma_n;\quad \sigma_n \sim t_\tau(0,1) & (5) \end{align}\] where \(t_\tau(0,1)\) is a standardized Student's t-distribution with \(\tau\) degrees of freedom (\(0 < \tau < 60\)). This choice accommodates fat-tailed returns, which are critical for modeling extreme events. The upper bound \(\tau < 60\) ensures the distribution remains sufficiently non-Gaussian while avoiding numerical instability (the t-distribution approaches the normal distribution as \(\tau \to \infty\)).

Regime-Dependent Volatility Shocks

To account for black swan events (e.g., the 2008 recession and 2020 pandemic), we introduce a state-dependent “amplitude” parameter modulating \(\mu_h\): \[\begin{align} \mu_h(t) = \begin{cases} \mu_h(t) + 0.8 \cdot \text{amplitude} & \text{if } t \in [262, 275] \quad \text{(2008 recession)}, \\ \mu_h(t) + 1.2 \cdot \text{amplitude} & \text{if } t \in [400, 410] \quad \text{(2020 pandemic)}, \\ \mu_h(t) & \text{otherwise.} \end{cases} \end{align}\] The multipliers (0.8 and 1.2) reflect the empirically observed difference in volatility magnitudes between these events (see Figure 2).

2.3.2 Simulation

T_breto_statenames <- c("H","G","Y_state")
T_breto_rp_names <- c("sigma_nu","mu_h","phi","sigma_eta", "tau", "amplitude")
T_breto_ivp_names <- c("G_0","H_0")
T_breto_paramnames <- c(T_breto_rp_names, T_breto_ivp_names)

T_breto_rproc1 <- "
  double beta,omega,nu,mu;
  mu = mu_h;
  if (t>=262 && t<=275){
      mu += 0.8 * amplitude;
    }
  else if (t>=400 && t<=410) {
      mu += 1.2 * amplitude;
  }
  omega = rnorm(0,sigma_eta * sqrt( 1- phi*phi ) * 
    sqrt(1-tanh(G)*tanh(G)));
  nu = rnorm(0, sigma_nu); // Change this
  G += nu; // Change This
  beta = Y_state * sigma_eta * sqrt( 1- phi*phi );
  H = mu*(1 - phi) + phi*H + beta * tanh( G ) 
    * exp(-H/2) + omega;
"

T_breto_rproc2.sim <- "
  double df = (nearbyint(tau) < 1) ? 1 : (nearbyint(tau) > 60) ? 60 : nearbyint(tau);
  Y_state = rt(df) * exp(H/2);  // t-distributed observation
 "

T_breto_rproc2.filt <- "
  Y_state = covaryt;
 "
T_breto_rproc.sim <- paste(T_breto_rproc1, T_breto_rproc2.sim)
T_breto_rproc.filt <- paste(T_breto_rproc1, T_breto_rproc2.filt)

T_breto_rinit <- "
  G = G_0;
  H = H_0;
  double df = (nearbyint(tau) < 1) ? 1 : (nearbyint(tau) > 60) ? 60 : nearbyint(tau);
  Y_state = rt(df) * exp(H/2);  // t-distributed observation
"

T_breto_rmeasure <- "
   y=Y_state;
"

T_breto_dmeasure <- "
  double df = (nearbyint(tau) < 1) ? 1 : (nearbyint(tau) > 60) ? 60 : nearbyint(tau);
  lik = dt(y/exp(H/2), df, give_log) - (H/2);  // Log-likelihood for t-distribution
"

T_breto_partrans <- parameter_trans(
  log=c("sigma_eta","sigma_nu"),
  logit="phi"
)

T_breto_filt <- pomp(data=data.frame(
    y=new_york_gasoline.ret.demeaned,time=1:length(new_york_gasoline.ret.demeaned)),
  statenames=T_breto_statenames,
  paramnames=T_breto_paramnames,
  times="time",
  t0=0,
  covar=covariate_table(
    time=0:length(new_york_gasoline.ret.demeaned),
    covaryt=c(0,new_york_gasoline.ret.demeaned),
    times="time"),
  rmeasure=Csnippet(T_breto_rmeasure),
  dmeasure=Csnippet(T_breto_dmeasure),
  rprocess=discrete_time(step.fun=Csnippet(T_breto_rproc.filt),
    delta.t=1),
  rinit=Csnippet(T_breto_rinit),
  partrans=T_breto_partrans
)

T_breto_params_test <- c(
  tau=5,
  sigma_nu = exp(-4.5),  
  mu_h = -5.5,
  amplitude = 2, 
  phi = expit(-0.5),     
  sigma_eta = exp(-10),
  G_0 = 0,
  H_0 = 0
)
  
T_breto_sim1.sim <- pomp(T_breto_filt, 
  statenames=T_breto_statenames,
  paramnames=T_breto_paramnames,
  rprocess=discrete_time(
    step.fun=Csnippet(T_breto_rproc.sim),delta.t=1)
)

T_breto_sim1.sim <- simulate(T_breto_sim1.sim,seed=1,params=T_breto_params_test)

T_breto_sim1.filt <- pomp(T_breto_sim1.sim, 
  covar=covariate_table(
    time=c(timezero(T_breto_sim1.sim),time(T_breto_sim1.sim)),
    covaryt=c(obs(T_breto_sim1.sim),NA),
    times="time"),
  statenames=T_breto_statenames,
  paramnames=T_breto_paramnames,
  rprocess=discrete_time(
    step.fun=Csnippet(T_breto_rproc.filt),delta.t=1)
)

We began our analysis with parameter simulation studies to initialize the model fitting process. The initial configuration: \[ \theta_0 = (\sigma_\nu,\mu_h,\phi,\sigma_\eta,G_0,H_0,\tau,\text{amplitude}) = (\exp(4.5), -5.5, \textrm{expit}(-0.5), \exp(-10),0,0,5,2) \] generated simulated trajectories that closely reproduced both the observed price dynamics and the amplified volatility during structural break periods (2008 recession at t=275 and 2020 pandemic at t=405; see Figure 8).

Using these initial values, we obtained a filtered log-likelihood of 457.797 (SE = 3.38e-6) for the simulated data. This initialization served as the basis for subsequent local optimization to refine parameter estimates.

2.3.3 Local Search

Figures 9 and 10 demonstrate satisfactory convergence in our local search, with log-likelihood values consistently increasing across all runs and remaining within a narrow range of a few log units. The key parametersrs (\(\mu_h, \phi, \tau\), and “amplitude”) exhibit stable convergence in most runs, though we observe a subset of runs with higher \(\sigma_\eta\) values that achieved marginally better likelihoods. The pair plot analysis reveals that our local search explored only a limited region of the parameter space, suggesting potential for further optimization. Notably, all runs now maintain ESS values above 345 in the final iteration, representing a substantial improvement over our previous model specification. To ensure robustness and more thoroughly investigate the parameter space, we conducted an additional global search with expanded parameter bounds.

2.3.4 Global Search

Our global search (Figures 11-12) achieved a maximum log-likelihood of 437.1, with values remaining relatively stable while showing consistent improvement across runs. The analysis reveals several key findings: First, all runs maintained effective sample size (ESS) values above 300 in the final iteration, demonstrating substantial improvement over previous model specifications. Second, while core parameters (\(\mu_h, \phi, \sigma_\eta\), and “amplitude”) converged stably in most runs, the pair plot indicates linear correlations between the log-likelihood and \((\log(\sigma_\nu), \mu_h, \phi, \sigma_\eta)\) This suggests the global search may not have fully explored the parameter space, as expanding the ranges for \(\phi\) and \(\sigma_\eta\) could potentially yield higher likelihood values. Although time constraints limited our ability to conduct more extensive searches, the current implementation nevertheless represents a significant advancement over our previous model configuration, both in terms of likelihood performance and sampling efficiency.

2.4 Modified Basic SV - No Leverage Model (POMP Framework)

2.4.1 Formulas

Our original hypothesis is that the leverage effect is more limited in gasoline prices compared to freely traded commodities or financial assets. To test this hypothesis, we fit a basic SV model without leverage. The basic SV model without leverage^[6][8][9] is defined as: \[\begin{align} Y_n &= \epsilon_n \exp(H_n / 2) \\ H_n &\sim N(\mu_h + \phi(H_{n-1} - \mu_h), \sigma) \\ H_1 &\sim N\left( \mu_h, \frac{\sigma}{\sqrt{1 - \phi^2}} \right) \end{align}\] where \(\epsilon_n \sim N(0,1)\). \(H_n\) is the latent parameter for the log volatility, \(\mu_h\) is the mean log volatility, and \(\phi\) is the persistence of the volatility term.

However, for the reasons specified above, we will model \(\epsilon_n \sim t_\tau(0,1)\) (\(0 < \tau < 60\)), and \(\mu_h\) as defined in (2.3.1).

2.4.2 Simulation

# Implement Basic Stochastic Model
T_basicSV_statenames <- "H"
T_basicSV_paramnames <- c("sigma","mu_h","phi", "tau", "amplitude")

T_basicSV_rproc <- "
  double mu;
  mu = mu_h;
  if (t>=262 && t<=275){
      mu += 0.8 * amplitude;
    }
  else if (t>=400 && t<=410) {
      mu += 1.2 * amplitude;
  }
  H = rnorm(mu + phi*(H - mu), sigma);
"

T_basicSV_rmeasure <- "
  double df = (nearbyint(tau) < 1) ? 1 : (nearbyint(tau) > 60) ? 60 : nearbyint(tau);
  y = rt(df) * exp(H/2);
"

T_basicSV_dmeasure <- "
  double df = (nearbyint(tau) < 1) ? 1 : (nearbyint(tau) > 60) ? 60 : nearbyint(tau);
  lik = dt(y/exp(H/2), df, give_log) - (H/2);
"

T_basicSV_rinit <- "
  H = rnorm(mu_h, sigma/(sqrt(1-phi*phi)));
"

T_basicSV_partrans <- parameter_trans(
  log=c("sigma"),
  logit="phi"
)

T_basicSV_filt <- pomp(data=data.frame(
    y=new_york_gasoline.ret.demeaned,time=1:length(new_york_gasoline.ret.demeaned)),
  statenames=T_basicSV_statenames,
  paramnames=T_basicSV_paramnames,
  times="time",
  t0=0,
  rmeasure=Csnippet(T_basicSV_rmeasure),
  dmeasure=Csnippet(T_basicSV_dmeasure),
  rprocess=discrete_time(step.fun=Csnippet(T_basicSV_rproc),
    delta.t=1),
  rinit=Csnippet(T_basicSV_rinit),
  partrans=T_basicSV_partrans
)

T_basicSV_params_test <- c(
  tau=5,
  sigma = exp(-4.5),  
  mu_h = -5.5,
  amplitude = 2,
  phi = expit(-0.1)
)

T_basicSV_sim <- simulate(T_basicSV_filt,seed=1,params=T_basicSV_params_test)

We initialized the estimation procedure with the following parameter configuration: \[ \theta_0 = (\sigma,\mu_h,\phi,\tau,\text{amplitude}) = (\exp(4.5), -5.5, \textrm{expit}(-0.1), 5,2) \] generated simulated trajectories that closely reproduced both the observed price dynamics and the amplified volatility during structural break periods.

Using these initial values, we obtained a filtered log-likelihood of 472.035 (SE = 6.99e-4) for the simulated data. We performed local search to refine parameter estimates.

2.4.3 Local Search

Figure 14 and 15 demonstrate satisfactory convergence in our local search, with log-likelihood values consistently increasing across all runs and remaining within a few log units of one another. However, one run achieved a notably lower likelihood, primarily due to suboptimal estimation of the persistence parameter \(\phi\).

The effective sample size (ESS) exhibits a significant drop around index 50, corresponding to the 1990 period. Despite this, the minimum ESS (280) represents an an acceptable value.

While key parameters (\(\mu_h, \phi, \sigma\)) show stable convergence, the pair plot reveals limitations in our local search, suggesting the explored parameter space may be too restricted. To address this, we conducted a subsequent global search to more thoroughly investigate the parameter space and ensure robustness in our estimates.

2.4.4 Global Search

All key parameters demonstrated stable convergence, with the log-likelihood showing consistent improvement across iterations. This robust convergence behavior provides strong support for the model's effectiveness in capturing the underlying price dynamics. The max log-likelihood value obtained was 434.8, with all runs maintaining effective sample size (ESS) values above 280 in the final iteration.

2.5 Model Comparison

To evaluate our hypothesis, we compared the performance of two competing models: the modified stochastic volatility (SV) model with leverage effects and its counterpart without leverage effects. Using the Akaike Information Criterion (AIC), defined as \(\text{AIC} = -2 \cdot l(\widehat{\theta}) + 2D\)^[10], where \(D\) represents the number of model parameters, we found the modified SV model without leverage effects to be statistically favored.

Model	LogLik	Parameters	AIC
Modified SV with Leverage	437.1	8	-858.2
Modified SV without Leverage	434.8	5	-859.6

This result aligns with our empirical observations from the global search of the leverage model (Section 2.3.4), where higher log-likelihood values consistently corresponded to smaller \(\log(\sigma_\nu)\) estimates (Figure 11). This relationship implies that the optimal value of \(\sigma_\nu\) approaches 0, causing \(G_n \approx G_{n-1}\) in equation (2), and consequently diminishing the leverage effect \(R_n\) in equation (1).

While these findings provide preliminary evidence supporting our hypothesis that gasoline prices exhibit more limited leverage effects than freely traded financial assets, we emphasize that more rigorous statistical testing would be required for definitive confirmation. The current analysis nevertheless offers valuable insights into the distinctive volatility dynamics of regulated commodity markets.

2.6 Model Comparison To The Benchmark T-GARCH

Lastly, we compared the modified SV model without leverage effects to a benchmark T-GARCH model. By incorporating the volatility as the latent state, we expect the SV model to outperform the T-GARCH model in terms of log-likelihood. Since we have been assuming a t-distribution, we fitted the T-GARCH model as a benchmark. The T-GARCH model^[6][8] with mean 0 is defined as:

\[\begin{align} Y_n &= \sigma_n \epsilon_n \quad \epsilon_n \sim t_\tau(0,1) \\ \sigma_n^2 &= \alpha_0 + \sum_{j=1}^p \alpha_j Y_{n-j}^2 + \sum_{k=1}^q \beta_k \sigma_{n-j}^2 \end{align}\]

The AIC table is shown below.

	q1	q2	q3	q4	q5	q6
p1	-861.402	-858.911	-865.018	-863.799	-862.378	-861.226
p2	-858.925	-857.911	-864.018	-862.799	-861.378	-860.226
p3	-857.921	-858.043	-863.018	-861.799	-860.378	-859.226
p4	-860.625	-859.625	-861.799	-860.799	-859.378	-858.269
p5	-859.661	-858.661	-860.378	-859.378	-858.378	-857.269
p6	-858.773	-857.773	-859.098	-858.098	-857.098	-856.269

## Fitting T-GARCH model with p = 1, q = 3

## Max log-likelihood: 435.509

The T-GARCH(3,1) model , selected via Akaike Information Criterion (AIC), achieved a higher log-likelihood (435.509) than our SV model (434.8), suggesting it may be better suited for capturing gasoline price volatility in this particular context. However, the relatively small margin of difference indicates that our SV framework still provides a reasonable approximation of the underlying dynamics, though with clear room for improvement.

Diagnostic checks (Figure 18) revealed important limitations in both approaches. Residual analysis showed the T-GARCH model's residuals deviated significantly from the assumed t-distribution, while autocorrelation functions revealed dependencies at lags 1 and 2. These patterns, combined with clear seasonal fluctuations identified through STL decomposition of the log returns (Figure 19), suggest that neither model adequately accounts for the full complexity of gasoline price dynamics. This represents a key limitation of our current analysis, as time constraints prevented us from incorporating seasonal components into our POMP framework.

3. Conclusion

Our investigation employed modified stochastic volatility (SV) models, both with and without leverage effects, to analyze gasoline price volatility dynamics. The model comparison, evaluated using the Akaike Information Criterion (AIC), revealed superior performance of the SV specification without leverage effects. This finding provides evidence that leverage effects in gasoline prices may be less pronounced than those observed in freely traded financial assets, potentially reflecting the unique characteristics of regulated energy markets. While these results support our initial hypothesis, we emphasize that additional robustness checks and formal hypothesis testing would be valuable to strengthen these conclusions. The current analysis nevertheless offers important insights into the distinct volatility patterns of gasoline prices compared to more conventional financial instruments.

4. Discussion

We made a serious mistake in specifying \(\mu_h\) the way we did in (2.3.1). Our original intension was to introduce a parameter to capture unexpected events. However, we ended up hard coding the period where the unexpected events happened. This is not a good practice since we basically introduce our biases into the model, which makes it incapable of adapting to new datasets. Given the time contraints, we propose the following approach to try out in the future. We suggest modeling event timing probabilistically rather than deterministically. Given the infrequent nature of such events, one could track elapsed time since the last simulated event using an accumulator variable, with event likelihood increasing over time. Perhaps this could be formalized through the use of a Weibull distribution.

Two additional limitations that are worth mentioning. First, we did not account for seasonal patterns in gasoline prices, which may have systematically influenced our volatility estimates. Incorporating seasonal components could improve model accuracy in future work. Second, our upper bound constraint of \(\tau < 60\) for the Student's t-distribution degrees of freedom was unnecessarily permissive. Since t-distributions with \(\tau > 30\) are practically indistinguishable from normal distributions, future implementations should enforce \(\tau < 30\) to preserve the intended fat-tailed properties while maintaining numerical stability. These refinements would strengthen the model's ability to capture extreme events without over-parameterization.

Prior work^[12] has applied Breto's stochastic volatility (SV) model with leverage effects to crude oil prices within the POMP framework, using a limited annual dataset spanning only 40 years. A key critique of that study was that the dataset might be insufficient for reliably estimating such a complex model, suggesting that a simpler specification could be more appropriate. This concern directly motivates our comparison between complex and simple models, as our own dataset — while larger — remains relatively short (~500 monthly observations).

The superior performance of a simpler model in our analysis might therefore reflect dataset limitations rather than invalidating our hypothesis about leverage effects. Indeed, a longer or higher-frequency (e.g., daily) series might reveal stronger evidence for leverage effects, as the additional data could better constrain the more complex model's parameters.

5. Contributions

Blinded.

6. Reference

[1] Natter, A. (2024, May 21). US puts its 1 million-barrel gasoline reserve up for sale. Bloomberg. https://www.bloomberg.com/news/articles/2024-05-21/us-puts-its-1-million-barrel-gasoline-reserve-up-for-sale

[2] Ghaddar, A., Lawler, A., & El Dahan, M. (2024, June 3). OPEC+ extends deep oil production cuts into 2025 | reuters. OPEC+ extends deep oil production cuts into 2025. https://www.reuters.com/business/energy/opec-seen-prolonging-cuts-2024-into-2025-two-sources-say-2024-06-02/

[3] https://www.eia.gov/dnav/pet/hist/LeafHandler.ashx?n=pet&s=eer_epmru_pf4_y35ny_dpg&f=m

[4] https://www.eia.gov/dnav/pet/TblDefs/pet_pri_spt_tbldef2.asp

[5] Mead, D., & Stiger, P. (2015, May). The 2014 plunge in import petroleum prices: What happened? U.S. BUREAU OF LABOR STATISTICS.

[6] Ionides, E. Lecture Notes for University of Michigan, STATS 531 Winter 2025. Modelling and Analysis of Time Series Data. Chapter 17.

[7] Bret'o C (2014). “On idiosyncratic stochasticity of financial leverage effects.” Statistics & Probability Letters, 91, 20-26. doi: 10.1016/j.spl.2014.04.003.

[8] Chan, J., & Grant, A. (2015). Modeling energy price dynamics: Garch Versus Stochastic volatility. SSRN Electronic Journal. https://doi.org/10.2139/ssrn.2616522

[9] https://mc-stan.org/docs/stan-users-guide/time-series.html#stochastic-volatility-models

[10] Ionides, E. Lecture Notes for University of Michigan, STATS 531 Winter 2025. Modelling and Analysis of Time Series Data. Chapter 5.

[11] https://www.rdocumentation.org/packages/lmenssp/versions/1.2/topics/qqplot.t

[12] https://ionides.github.io/531w22/final_project/project18/blinded.html

[13] UM ChatGPT was used to polish the sentences and correct grammars.