Updating model_gbm to accurately estimate mu

[maqifrnswa]

Awesome write up, thanks - I learned a lot!

I found why the MCMC result couldn't find the right mean if you set it to a non-zero value. The current model_gbm is

model_gbm = pm.Model()

with model_gbm:
    prior_mean = pm.Normal('mean', mu=prior_mean_mu, sd=prior_mean_sigma)
    prior_vol = pm.HalfNormal('volatility', sd=prior_vol_sigma)
    
    likelihood = pm.Normal(
        'daily_returns', mu=prior_mean, sd=prior_vol, observed=returns_geo)

When the returns_geo is log (St/S0), the mean of the likelihood's mean is prior_mean - 0.5 * prior_vol **2. Right now, the code finds the likelihood distribution mean and in post processing tries to find the value of the true mean, but that results in the fitting not accounting for variance in the return. If you do the following instead:

model_gbm=pm.Model()
with model_gbm:
    prior_mean = pm.Normal('mean', mu=prior_mean_mu, sigma=prior_mean_sigma)
    prior_vol = pm.HalfNormal('volatility', sd=prior_vol_sigma)
    drift = pm.Deterministic("log drift", prior_mean-0.5 *prior_vol**2 )
    likelihood  = pm.Normal("daily returns", mu=drift, sigma=prior_vol, observed=returns_geo_full)

that deterministic random variable (which might be a newer feature in pymc) will account for the proper stochastic mixing of variance and mu, and it will accurately estimate mu from the data.