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).
Vasicek drift rate: Parameters $\hat{\pi}$ (long-term mean), $\hat{\kappa}$ (mean reversion), $\hat{\omega}$ (volatility).
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).
Vasicek drift rate: $\hat{\pi}$ = 0.030839 (0.003630), $\hat{\kappa}$ = 9.128286 (0.633592), $\hat{\omega}$ = 0.332236 (0.009359).
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!
