Ожидаемая цена акций не одинакова, и зависит от волатильности

Сделал повторный фиттинг на специально подготовленных исторических данных, исключив по возможности предвзятости и неравномерности, 250 акций, 50 лет истории каждая (хотя, bias — ошибка выжившего остается). Получается сильная зависимость от волатильности.

Например, ожидаемая годовая прибыль получается для INTC х1.18, KO х1.06.

Модель для фиттинга E[R_365] = a + b*r_rf_365 + c*sqrt(var_365) оптимизация наименьших квадратов, все переменные в оригинальном, линейном пространстве, кроме волатильности она посчитана в лог маштабе, как var_365 = 0.8*252*EWMA[log r_daily].

Доверительные интервалы хорошие

a = 0.545548  (95% CI: 0.486466, 0.604631)
b = 0.423923  (95% CI: 0.367767, 0.480078)
c = 0.502396  (95% CI: 0.489426, 0.515367)

Проблема: насколько я понимаю это противоречит классической теории которая говорит что волатильность не влияет на ожидаемую цену. И это противоречит интуиции, мне кажется годовая прибыль х1.18 для Интела и тому подобных волатильных акций слишком высока, я интуитивно оценил бы ее ниже, как 1.1. Но, ошибки вроде нет… я также сделал другой фиттинг, с оптимизацией Likelihood в лог пространстве, он дал похожие результаты.

Если есть желание проверить:

Подготовленные данные с уже расчитанными волатильностями и код

# Fitting Risk Premium from Historical Data
#
#   E[R_T] = a + b*r_rf_T + c*sqrt(var_T)

import numpy as np
import statsmodels.api as sm
import numpy as np
import pandas as pd

period = 365

def load_data():
  df = pd.read_csv('/storage/data/alien/experiments/stock_prices/fitting_data.tsv', sep='\t')
  df['ar_t'] = np.exp(df['r_t'])
  return df

def fit_risk_premium(df, period_d):
  df = df[df['period'] == period_d].copy()

  ar_rf_1y_t0 = df['ar_rf_1y_t0'].values
  ar_rf_T = np.exp((period_d / 365) * np.log(ar_rf_1y_t0))
  y = df[f'ar_t'].values

  # 0.69 - trading days multiplier, 0.8 - some constant, optional
  var_T = df['mvar_d_t0'].values * period_d * 0.69 * 0.8

  X = np.column_stack([
    ar_rf_T,
    np.sqrt(var_T),
  ])

  X = sm.add_constant(X)
  model = sm.OLS(y, X, missing='drop').fit()

  a = model.params[0]
  b = model.params[1]
  c = model.params[2]

  def print_info():
    ci = model.conf_int()
    mse = np.mean((model.fittedvalues - model.model.endog) ** 2)
    print(f"Fitting model: E[AR_{period_d}] = a + b*ar_rf_T + c*sqrt(var_T) + d*var_T")
    print(f"{'R²':>10} = {model.rsquared:.6f}")
    print(f"{'MSE':>10} = {mse:.6f}\n")
    print(f"{'a':>10} = {a:.6f}  (95% CI: {ci[0][0]:.6f}, {ci[0][1]:.6f})")
    print(f"{'b':>10} = {b:.6f}  (95% CI: {ci[1][0]:.6f}, {ci[1][1]:.6f})")
    print(f"{'c':>10} = {c:.6f}  (95% CI: {ci[2][0]:.6f}, {ci[2][1]:.6f})")

  print_info()

  def risk_premium(ar_rf_1y_t0_input, var_input):
    ar_rf_input = np.exp((period_d / 365) * np.log(np.array(ar_rf_1y_t0_input)))
    var_input = np.array(var_input) * period_d * 0.69 * 0.8
    return (
      a + b * ar_rf_input + c * np.sqrt(var_input)
    )

  return (a, b, c), risk_premium

def run():
  period = 365
  df = load_data()

  ks, risk_premium = fit_risk_premium(df, period)
  print('INTC', risk_premium(1.0505, 0.0007))
  print('KO', risk_premium(1.0505, 0.0001))


run()

И, чистые исходные данные, цены акций и безрисковой ставки если будет желание использовать другую меру волатильности и посчитать ее самостоятельно. Я использовал EWMA(span=365/3)