估計不同採樣率下幾何布朗運動的波動性
當以不同的採樣頻率重新採樣數據時,我難以估計波動率(= 對數回報的標準偏差)。
問題
我使用幾何布朗運動生成了時間序列數據。原始時間序列數據以1 小時間隔生成半年:
$$ \begin{align} \mathrm{days} & = 365 / 2 = 182.5,, \ n & =182.5 \cdot 24 = 4380,, \ \Delta t & = \mathrm{days} / 365 / n=0.000114155,. \end{align} $$
假設:
- 全年交易連續,交易日數為365天。
- $ t=1 $ 表示 1 年(=365 個交易日)。
原始 GBM 時間序列 $ {x_1 , \ldots , x_n} $ 使用以下參數生成:
$$ \begin{align} \mu & = 0.5,, \ \sigma_\mathrm{1h} & = 0.8,, \ x_0 & = 1000,, \ \Delta t & = 0.000114155,. \end{align} $$
使用此參數,日波動率和年波動率分別為:
$$ \begin{align} \sigma_\mathrm{daily} & = \sigma_\mathrm{1h} \sqrt{24} = 3.92,,\ \sigma_\mathrm{annual} & = \sigma_\mathrm{daily} \sqrt{365} = 74.88,. \ \end{align} $$
參數估計
現在我使用以下公式估算波動率:
$$ \begin{align} r_i & = \log{\frac{x_i}{x_{i-1}}} \qquad \textrm{log returns}\ \hat{\mu} & = \frac{1}{n}\cdot\sum_{i=1}^n{r_i} \qquad \textrm{mean return}\ S_\mu^2 & = \frac{1}{n-1}\sum_{i=1}^n{(r_i-\hat{\mu})^2} \qquad \textrm{sample variance of returns}\ \hat{\sigma} & = \frac{S_\mu}{\sqrt{\Delta t}} \qquad \textrm{estimated volatility} \end{align} $$
將其應用於原始時間序列(1h 間隔)數據,我得到正確的結果:
$$
\begin{align} \hat{\sigma}\mathrm{1h} & = 0.7970 \ \hat{\sigma}\mathrm{daily} & = 3.871 \ \hat{\sigma}_\mathrm{annual} & = 73.959 \ \end{align} $$從重新採樣的數據估計參數
但是,當我將時間序列數據重新採樣為1 天間隔(或任何其他間隔)時,估計的 sigma $ \hat\sigma $ 保持不變(大約 0.8),但由於頻率不同, $ \Delta t $ 是不同的並且 $ \hat{\sigma}\mathrm{daily} $ 和 $ \hat{\sigma}\mathrm{annual} $ 是不正確的。
將時間序列重新採樣到1 天箱,給出
$$
\begin{align} n & = 183 \ \Delta t & = \textrm{days} / 365 / n = 0.002725 \ \hat{\sigma}\mathrm{1d} & = 0.773 \ \hat{\sigma}\mathrm{daily} & = 0.775 \ \hat{\sigma}_\mathrm{annual} & = 14.811 \ \end{align} $$我很確定問題出在波動率的某個地方,但我無法弄清楚這裡到底出了什麼問題。
重新採樣的(1 天箱)數據在每個數據點之間確實有不同的時間間隔,所以 $ \Delta t $ 必須不同。
用 Python 實現
import numpy as np import matplotlib.pyplot as plt np.random.seed(1018) mu = 0.5 sigma = 0.8 x0 = 1000 start = pd.to_datetime("2021-01-01") end = pd.to_datetime("2021-07-02 11:59.9999") freq = "1H" date_index=pd.date_range(start, end, freq=freq) n = len(date_index) dur = end-start dur_days = (dur.days + dur.seconds / (3600 * 24)) dt = dur_days / 365 / n print("n:", n) print("dur_days: %.1f" % dur_days) print("dt: %.8f" % dt) x = np.exp((mu - sigma ** 2 / 2) * dt + sigma * np.random.normal(0, np.sqrt(dt), size=n).T) x = x0 * x.cumprod(axis=0) df0 = pd.DataFrame({"values": x}, index=date_index) df = df0.copy() real_daily_vola = sigma * np.sqrt(1 / dt / 365) real_ann_vola = real_daily_vola * np.sqrt(365) print("sample vola (=sigma): %.2f for freq %s" % (sigma, freq)) print("real daily vola: %.2f" % real_daily_vola) print("real ann vola: %.2f" % real_ann_vola) def estimate_vola(s): n = len(s) start = s.index[0] end = s.index[-1] dur = end - start dur_days = dur.days + dur.seconds / (3600 * 24) dt = dur_days / 365 / n r = np.log(s / s.shift(1)) mu = np.nanmean(r) s2 = 1 / (n-1) * ((r - mu)**2).sum() est_sigma = np.sqrt(s2 / dt) daily_vola = est_sigma * np.sqrt(1 / dt / 365) ann_vola = daily_vola * np.sqrt(365) return est_sigma, daily_vola, ann_vola print("Estimate original vola: sigma=%.2f daily_vola=%.2f ann_vola=%.2f" % estimate_vola(df0["values"])) # Re-sample to 1d interval df = pd.DataFrame({"open": x, "high": x, "low": x, "close": x}, index=date_index) df_1d = df.resample('1D').agg({'open': 'first', 'high': 'max', 'low': 'min', 'close': 'last'}) print("Re-sampled to 1d bins: sigma=%.2f daily_vola=%.2f ann_vola=%.2f" % estimate_vola(df_1d["close"])) # Re-sample to 2d interval df = pd.DataFrame({"open": x, "high": x, "low": x, "close": x}, index=date_index) df_2d = df.resample('2D').agg({'open': 'first', 'high': 'max', 'low': 'min', 'close': 'last'}) print("Re-sampled to 2d bins: sigma=%.2f daily_vola=%.2f ann_vola=%.2f" % estimate_vola(df_2d["close"]))
輸出
n: 4380 dur_days: 182.5 dt: 0.00011416 sample vola (=sigma): 0.80 for freq 1H real daily vola: 3.92 real ann vola: 74.88 Estimate original vola: sigma=0.79 daily_vola=3.87 ann_vola=73.96 Re-sampled to 1d bins: sigma=0.77 daily_vola=0.78 ann_vola=14.81 Re-sampled to 2d bins: sigma=0.79 daily_vola=0.56 ann_vola=10.70
筆記
1 / dt / 365
線性比例因子daily_vola = est_sigma * np.sqrt(1 / dt / 365)
轉換 $ \hat\sigma $ 我認為進入日常 vola 是正確的。這是:
重新採樣頻率 1 / dt / 365
一維 1 二維 0.5
您生成路徑的方式,您期望年化vol 為 0.8(因為您乘以 $ \sqrt{dt} $ , 和 $ dt $ 以年表示。
在實現的計算中,您可以使用原始數據將其取回
np.std(np.diff(np.log(df0['values']),1))*np.sqrt(365*24)
或在重新採樣的每日數據上分別使用
np.std(np.diff(np.log(df_1d["close"]),1))*np.sqrt(365) np.std(np.diff(np.log(df_2d["close"]),1))*np.sqrt(365/2)
正如預期的那樣,所有 3 個都給出了接近但不完全是 0.8 的結果。