Replicating Maximum-Likelihood Inflation Model Estimates

# Replicating Parameter Estimates from Fergusson’s 2020 Inflation Forecasting Models


# Replicating Parameter Estimates from Fergusson’s 2020 Inflation Forecasting Models












I’m working on replicating the parameter estimates for the four models in Table 1 of the paper “Forecasting Inflation Using Univariate Continuous-Time Stochastic Models” by Kevin Fergusson (2020), available to download at ResearchGate. The models, estimated using maximum likelihood estimation (MLE) on US CPI data (from FRED, January 1913 to July 2018), are:

- Log-normal: Parameters $\hat{\pi}$ (drift) and $\hat{\sigma}$ (volatility).

- Log-normal with local drift: Parameter $\hat{\sigma}$ (volatility, with drift inferred via filter).



- Log-normal with Vasicek drift rate: Parameters $\hat{\pi}$, $\hat{\kappa}$, $\hat{\omega}$, and $\hat{\sigma}$.

Table 1 reports the following estimates (with standard errors):

- Log-normal: $\hat{\pi}$ = 0.031053 (0.002186), $\hat{\sigma}$ = 0.022439 (0.000446).

- Log-normal with local drift: $\hat{\sigma}$ = 0.021003 (0.000418).



- Log-normal with Vasicek drift rate: $\hat{\pi}$ = 0.031122 (0.006782), $\hat{\kappa}$ = 1.745702 (0.195865), $\hat{\omega}$ = 0.121449 (0.002592), $\hat{\sigma}$ = 0.019800 (0.000424).

I’m using Python and have sourced the CPI data from FRED, but I’m struggling with replicating the parameter estimates, particularly for the last 3 models. I was able to replicate the Log-normal model parameter estimates almost perfectly, but for the drift models my challenge is coding the filter (stated in section 4.2) used to get the inflation rates ${\pi}_{t_i}$.

I have python code for the Log-normal model and Vasicek drift rate model, but I think showing both of those in an excel is easier. See the Google Drive Excel for those calculations. The Log-normal model estimates are almost perfect while the Vasicek drift rate model estimates have small differences, but since we are using historical data there should be no differences.

Has anyone replicated these parameter estimates or implemented similar continuous-time stochastic models? Could you share Python (or R) code for the MLE process?

For instance, for the Log-normal with local drift model, I tried recreating the filter from the paper:

> …where we infer the unobserved rates by applying a filter f to the observed noisy logarithms of the index log($I_{t_j}$), j ∈ {i−11, i−10, …, i}, and computing ${\pi}_{t_i}$ as the first-order derivative of f at time ${t_i}$. The filter f is obtained by an exponentially weighted regression of log($I_{t_j}$) on the basis terms of the quadratic $𝑓(t_𝑗) = a + bt_𝑗 + ct_{𝑗}^2$ 𝑗 for j ∈{i − 11, i − 10, …, i}, from which the drift rate is given by $\pi_{t_i} = \frac{df(t)}{dt} \big|_{t=t_i} = b + 2c t_i$.

We use the unobserved rates from the filter to calculate the parameter: $$ \hat{\sigma}^2 = \frac{2}{\Delta} \left[ -1 + \sqrt{1 + \frac{1}{n} \sum_{i=1}^{n} \left( \log I_t - \log I_{t-1} - \pi_{t-1} \Delta \right)^2} \right] $$ Below is my attempt for the Log-normal with local drift model:

```
import pandas as pd
import numpy as np

# Load data from Excel file
excel_file = '2020 - Fergusson - Data.xlsx'
df = pd.read_excel(excel_file, sheet_name='Data')
# Convert observation_date to datetime and sort
df['observation_date'] = pd.to_datetime(df['observation_date'])
df = df.sort_values('observation_date')

# Compute log CPI and monthly differences
df['log_CPI'] = np.log(df['CPIAUCNS'])
df['diff_log_CPI'] = df['log_CPI'].diff()

# Set time variable in years
m = len(df)
df['t'] = (np.arange(m) / 12)

# Parameters
Delta = 1/12  # Monthly time step
lambda_ = 0.25  # Decay parameter for exponential weighting

# Estimate sigma_hat with exponentially weighted regression
sum_A = 0
count = 0
for i in range(12, m):
    # Use a 12-month window of past data
    window = df.iloc[i-12:i]
    t_window = window['t'].values
    log_I_window = window['log_CPI'].values
    # Exponential weights: w_j = exp(-lambda * j)
    weights = np.exp(-lambda_ * np.arange(12)[::-1])
    # Fit quadratic polynomial with weights
    coeffs = np.polyfit(t_window, log_I_window, 2, w=weights)
    c, b, a = coeffs
    # Drift π_{t-1} is the derivative at t[i-1]: b + 2*c*t[i-1]
    t_im1 = df['t'].iloc[i-1]
    pi_t_im1 = b + 2 * c * t_im1
    # Residual
    r_i = df['diff_log_CPI'].iloc[i]
    sum_A += (r_i - pi_t_im1 * Delta)**2
    count += 1

# Compute sigma_hat using the formula
A = sum_A / (count)
sigma2_hat = (2 / Delta) * (-1 + np.sqrt(1 + A))
sigma_hat = np.sqrt(sigma2_hat)

print(f"Estimated volatility (σ̂): {sigma_hat:.6f}")
#Output: Estimated volatility (σ̂): 0.021543
#Target is 0.021003
```

I haven't attempted the Log-normal with Vasicek drift rate since it wouldn't make sense until I get the other models correct.

Any help for any of the parameter estimates for any of the models would be much appreciated! Thanks for any insights or code samples!