风管模型第四次作业

ErnestDong · Jun 7, 2023 · b429593 · b429593
1 parent 301acab
commit b429593
Show file tree

Hide file tree

Showing 6 changed files with 1,597 additions and 1,236 deletions.
diff --git a/main.ipynb b/main.ipynb
diff --git a/main.py b/main.py
@@ -1,215 +1,185 @@
 # %% [markdown]
-# # 风管模型第三次作业
+# # 风管模型第四次作业
 #
+# 股票市场的走势是一个国家经济状况的“晴雨表”。长久以来，此领域一直是众多学者的研究重点。目前，预测股票价格走势的方法多种多样，但是均存在对股票价格的波动拟合效果较差、预测精度有限等问题。时间序列模型具有应用范围广、限制要求少、短期预测准确率高等优点，借此，它成为了金融预测领域较受欢迎的预测模型之一。
 #
-
-# %% [markdown]
-# ## 单个公司的 KMV 模型
-#
-# ### 假设
-#
-# - 市场的无风险利率为3%
-# - 公司的违约实施点为企业1年以下短期债务的价值加上未清偿长期债务账面价值的一半
-# - 贷款期限一年
-# - 当公司发生违约时，违约回收率为0，即全部贷款均不能回收。
-# - 贷款利率为包括流动性溢价和信用溢价，分别用无风险利率和预期违约率$\times$100表示，即贷款利率 = 无风险利率 + 预期违约率 $\times$ 100
-#
-# ### 计算股票价格波动率
-#
-# 我们选取了三家上市公司，分别为比亚迪，三一重工，宁德时代，比亚迪为汽车厂商，三一重工为国有大型装备制造厂商，宁德时代为我国新能源电池领头企业，三家厂商都具备较好的现金流以及信用背书，因此是我们选取的贷款对象。
+# 时间序列分析是一种主要用于描述和探究某一现象随时间的发展如何变化的分析方法。时间序列是同一种现象在不同时间上的相继观察值排列而成的一组数字序列。人们运用时间序列的历史数据，揭示现象随时间变化的规律并将这种规律延伸到未来，从而对该现象的未来做出预测。
 #
-# 根据三家公司的历史股价数据，我们计算出三家厂商的股票波动率。
+# 本文运用时间序列分析方法，使用LN、AR、ARCH、GARCH、RSLN-2五个模型对样本数据进行拟合，计算对数似然函数的最大值以及AIC和SBC，比较不同模型的拟合效果，并对沪深300指数未来的走势进行预判。
 
 # %%
-import glob
-
+import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 from arch import arch_model
-from scipy import optimize, stats
+from statsmodels.tsa.arima.model import ARIMA
+from statsmodels.tsa.regime_switching.markov_regression import MarkovRegression
 
-companies = glob.glob("data/trade/*.xlsx")
+# print("importing done")
+plt.rcParams['font.sans-serif']=['SimHei']
+plt.rcParams['axes.unicode_minus']=False
 result = {}
-rate = {}
-for company in companies:
- df = pd.read_excel(company, index_col=0, skiprows=1).dropna()
- company = company.split("/")[-1].split("[")[0]
- result[company] = {
- "V": df["总市值"].apply(lambda x: float(x.replace(",", ""))) * 10000
- }
- rate[company] = df["涨跌幅"] / 100
-rate = pd.DataFrame(rate).sort_index()
-rate
+df = (
+ pd.read_excel("./lib/创业板指.xlsx", index_col=2, skipfooter=2)
+ .sort_index()
+ .drop(columns=["代码", "名称"])
+)
+df = df.apply(
+ lambda x: x.str.replace(",", "").astype(float) if x.dtype == "object" else x
+)
+df.index = pd.to_datetime(df.index)
+df["weeknum"] = df.index.map(lambda x: str(x.year) + "-" + str(x.isocalendar().week))
+df = df.sort_index()
+df["log(收益率)"] = (df["收盘价(元)"] / df["收盘价(元)"].shift(1)).dropna().apply(np.log)
+df[["log(收益率)"]].plot(kind="kde")
 
 # %% [markdown]
-# ### 利用 garch 模型计算相关系数与协方差阵
+# ## 第一题
+#
+# ### LN 模型
+#
+# 我们依据<cite data-cite="wang2020corrected">Wang(2020)</cite> 提出的算法，对 Lognormal 模型进行最大似然估计，即
 #
-# 假设三只股票的收益率服从GARCH (1,1)，然后算出得到标的资产组合的协方差矩阵和相关阵。
+# $$
+# \log(L) = -\sum_{i=1}^n\log x_i-n\log\sigma-\frac{n}{2}\log 2\pi - \sum_{i=11^n}\frac{(\log x_i-\mu)^2}{2\sigma^2}
+# $$
+#
+# 采用样本均值 $\hat\mu$ 和样本方差 $\hat\sigma^2$ 代入计算 $AIC$、$SBC$，得到的结果如下表所示
 
 # %%
-cov = {}
-var = {}
-r2 = {}
-rate = rate.dropna()
-for company in companies:
- company = company.split("/")[-1].split("[")[0]
- r2[company] = {}
- cross = rate[company] * rate[company] * 1000
- arch = arch_model(cross)
- arch_param = arch.fit(update_freq=0)
- var[company] = (arch_param.forecast(reindex=True).variance.values[-1][0]) / 1000
- parameters = arch_param.params[1:]
- cross = rate[[i for i in rate.columns if i != company]].prod(axis=1)
- covlist = [cross[0]]*100000
- print(parameters)
- for j in range(1, len(cross)):
- covlist.append(parameters[0] + parameters[1] * cross[j] + parameters[2] * covlist[-1])
- cov[company] = covlist[-1]/100000
-covar = {i: {} for i in cov}
-for company in cov:
- r2[company][company] = 1
- covar[company][company] = var[company]
- others = tuple(i for i in cov if i != company)
- r2[others[0]][others[1]] = cov[company] / (
- var[others[0]] ** 0.5 * var[others[1]] ** 0.5
- )
- r2[others[1]][others[0]] = r2[others[0]][others[1]]
- covar[others[0]][others[1]] = cov[company] ** 0.5
- covar[others[1]][others[0]] = cov[company] ** 0.5
-corr = pd.DataFrame(r2)
-cov = pd.DataFrame(covar)
-corr
+# df["weeknum"] = df.index.map(lambda x: str(x.year) + "-" + str(x.isocalendar().week))
+# df = df.drop_duplicates(subset="weeknum", keep="last")
+log = (df["收盘价(元)"] / df["收盘价(元)"].shift(1)).dropna().apply(np.log)
+mu = log.mean()
+sigma = log.std()
+likelyhood_ln = log.apply(
+ lambda x: -np.log(sigma) - np.log(2 * np.pi) / 2 - (x - mu) ** 2 / (2 * sigma**2)
+).sum()
+k, logl = 2, likelyhood_ln
+result["LN"] = {
+ "k": k,
+ "log(L)": logl,
+ "AIC": logl - k,
+ "SBC": logl - k * np.log(log.shape[0]) / 2,
+}
+pd.DataFrame(result).T
 
 # %% [markdown]
-# ### 计算违约点
+# ### AR(1) 模型
+#
+# 自回归模型（AR）是最简单的时间序列处理方法，对LogNormal的改进在于考虑了变量之间的相关性，用过去随机变量的情况来描述未来该随机变量的取值。在自回归模型假设下，随机变量的取值总会向长期均值靠拢。我们选择一阶自回归模型对样本数据进行拟合，函数表达式为
+# $$Y_t=\mu+\alpha(Y_{t-1}-\mu)+\sigma\varepsilon_t$$
+# 其中 $\alpha\in(-1,1)$ ，$\varepsilon_t$ 为相互独立且符合标准正态分布的白噪声过程。
 
 # %%
-companies = glob.glob("data/balance/*.xlsx")
-for company in companies:
- df = (
- pd.read_excel(company, index_col=0, skipfooter=4, sheet_name="file")
- .T.rename(columns=lambda x: x.strip())
- .head(1)
- )
- company = company.split("/")[-1].split("(")[0]
- result[company]["F"] = (df["流动负债合计"] + 0.5 * df["非流动负债合计"]).mean()
-{company:result[company]["F"] for company in result}
+model = ARIMA(log.values, order=(1, 0, 0))
+AR1 = model.fit()
+k, logl = 3, AR1.llf
+result["AR(1)"] = {
+ "k":k,
+ "log(L)": logl,
+ "AIC": logl - k,
+ "SBC": logl - k * np.log(log.shape[0]) / 2,
+}
+print(AR1.summary())
+pd.DataFrame(result).T
 
 # %% [markdown]
-# ### 计算预期违约概率 EDF
+# 自回归方法的优点是所需资料不多，可用自身变数数列来进行预测。但是这种方法受到一定的限制：必须具有自相关，自相关系数是关键。如果自相关系数(R)小于0.5，则不宜采用，否则预测结果极不准确。自回归只能适用于预测与自身前期相关的经济现象，即受自身历史因素影响较大的经济现象，如矿的开采量，各种自然资源产量等；对于受社会因素影响较大的经济现象，不宜采用自回归，而应改采可纳入其他变数的向量自回归模型。
+#
+# ### ARCH 模型
 #
-# 利用KMV模型，我们计算了三家公司各自的预期违约概率，结果如下：
+# ARCH模型即自回归条件异方差模型，是一种用来处理时间序列的模型。在股票中，ARCH可以用来预测股票的波动率，从而控制风险。粗略地说，该模型将当前一切可利用信息作为条件，并采用某种自回归形式来刻画方差的变异，对于一个时间序列而言，在不同时刻可利用的信息不同，而相应的条件方差也不同，利用ARCH模型，可以刻画出随时间而变异的条件方差。ARCH模型解决了LogNormal和AR被诟病的波动率恒定问题。这里我们同样使用一阶条件异方差模型对样本数据进行拟合，函数表达式为 $$Y_t=\mu+\sigma_t \varepsilon_t$$
+# 其中 $\sigma_t^2=a_0+a_1(Y_(t-1)-\mu)^2$。
 
 # %%
-def kmv(r, sigma_e, t, equity, debt):
- def option(w):
- x, sigma_a = w
- N_d1 = stats.norm.cdf(
- (np.log(abs(x) * equity / debt) + (r + 0.5 * sigma_a**2) * t)
- / (sigma_a * np.sqrt(t))
- )
- N_d2 = stats.norm.cdf(
- (np.log(abs(x) * equity / debt) + (r - 0.5 * sigma_a**2) * t)
- / (sigma_a * np.sqrt(t))
- )
- e1 = equity - (x * equity * N_d1 - debt * N_d2 * np.exp(-r * t))
- e2 = sigma_e - sigma_a * N_d1 * x
- return [e1, e2]
-
- assets, sigma_a = optimize.fsolve(option, [1, 0.1], maxfev=1000000000)
- DD = (assets * equity - debt) / (assets * equity * sigma_a)
- EDF = stats.norm.cdf(-DD)
- return EDF, DD
-
-
-edf = {}
-for company in result:
- kmvmodel = kmv(
- 0.03,
- result[company]["V"].pct_change().var() * 252,
- 1,
- result[company]["V"].values[0],
- result[company]["F"],
- )
- edf[company] = {"edf": kmvmodel[0], "dd": -kmvmodel[1]}
-edf = pd.DataFrame(edf)
-edf
+arch = arch_model(log, p=1, q=0, rescale=False).fit(update_freq=0)
+k, logl = 3, arch.loglikelihood
+print(arch.summary())
+result["ARCH"] = {
+ "k":k,
+ "log(L)": logl,
+ "AIC": logl - k,
+ "SBC": logl - k * np.log(log.shape[0]) / 2,
+}
+pd.DataFrame(result).T
 
 # %% [markdown]
-# ## 组合的违约概率
+# ### GARCH 模型
 #
-# 依据<cite data-cite="任宇航2006信用风险组合管理模型中的相关性问题研究述评">任宇航(2006)</cite>提出的算法，$DD$ 服从联合多元标准正态分布，且相关系数等于资产的相关系数，因而以下 code cell 分别为全部违约、两个公司违约、三个公司违约的概率
+# GARCH模型，即广义ARCH模型，它是对ARCH模型的进一步扩展。它同时考虑了序列相关和条件方差，对误差的方差进行了进一步建模，具有更强的灵活性。本文选择 $GARCH(1,1)$ 模型对历史数据进行拟合。
+#
+# 根据拟合结果可知，对数似然函数值为 7557.05，AIC为 7553.0492，SBC为 7541.0919。
 
 # %%
-cov = rate.cov()
-alldefault = stats.multivariate_normal.cdf(x=edf.T["dd"], cov=cov)
-alldefault
+garch = arch_model(log, p=1, q=1, rescale=False).fit(update_freq=0)
+k, logl = 4, garch.loglikelihood
+print(garch.summary())
+result["GARCH"] = {
+ "k":k,
+ "log(L)": logl,
+ "AIC": logl - k,
+ "SBC": logl - k * np.log(log.shape[0]) / 2,
+}
+pd.DataFrame(result).T
 
-# %%
-twodefault = {}
-for company in edf:
- others = [i for i in edf if i != company]
- others_edf = edf[others]
- others_cov = cov[others].loc[others]
- twodefault[tuple(others)] = (
- stats.multivariate_normal.cdf(x=others_edf.T["dd"], cov=others_cov.values)
- - alldefault
- )
-twodefault
+# %% [markdown]
+# ### RSLN-2
+#
+# RSLN 模型本质是马尔可夫状态转移模型，因而我们参考[statsmodels 文档](https://www.statsmodels.org/dev/generated/statsmodels.tsa.regime_switching.markov_regression.MarkovRegression.html)，令两 regime 分布异方差，回归结果如下表，显示异方差的确存在(不同状态下的方差影响显著)。
+#
+# 同样，我们提取似然函数对数，并计算 $AIC$、$SBC$
 
 # %%
-onedefault = {}
-for company in edf:
- others = [i for i in edf if i != company]
- onedefault[company] = (
- edf.T["edf"].to_dict()[company] - alldefault - twodefault[tuple(others)]
- )
-onedefault
+rsln = MarkovRegression(log.values, 2, switching_variance=True).fit() # 异方差
+k, logl = 6, rsln.llf
+print(rsln.summary())
+result["RSLN-2"] = {
+ "k": k,
+ "log(L)": logl,
+ "AIC": logl - k,
+ "SBC": logl - k * np.log(log.shape[0]) / 2,
+}
+pd.DataFrame(result).T
 
 # %% [markdown]
-# ## 贷款比例与信用风险
+# ### 综合比较
 #
-# 利用以上的结果，我们以贷款给这三家公司的期望收益为标准，选择期望收益最大化的贷款分配比例。由于违约概率都很小，贷款的 $VaR$ 指导意义较弱，因而我们选择最大化期望收益。最终的结果是100%宁德时代，以获得最高的期望收益。
+# 不同似然函数对比采用 $\chi^2$ 检验，以 lognormal 模型为基准进行 LRT 检验，结果如下表显示。综合似然函数最大值、$AIC$、$SBC$ 与 $LRT$ 检验值，GARCH 模型的似然函数、$AIC$、$SBC$ 最大，且相较于 lognormal 模型效果提升显著。
 
 # %%
-zerodefault = 1 - alldefault - sum(twodefault.values()) - sum(onedefault.values())
-return_rate = (1.03+edf.T["edf"] * 10**3).to_dict()
-def expected_return(case):
- case = {'比亚迪': case[0], '三一重工':case[1], '宁德时代': case[2]}
- final = (
- alldefault * 0
- + sum(
- [
- twodefault[i] * return_rate[j] * case[j]
- for i in twodefault
- for j in onedefault
- if j not in i
- ]
- )
- + sum(
- [
- onedefault[i] * return_rate[j] * case[j]
- for i in onedefault
- for j in onedefault
- if j != i
- ]
- )
- + zerodefault * sum([return_rate[i] * case[i] for i in onedefault])
- )
- return -final
-
-
-case = optimize.minimize(
- expected_return,
- (0,0,1),
- constraints=({"type": "eq", "fun": lambda x: sum(x) - 1}),
- bounds=[(0, 1)] * 3,
-).x
-{'比亚迪': case[0], '三一重工':case[1], '宁德时代': case[2]}
+# basemodel = "LN"
+# for model in result:
+# if model == basemodel:
+# continue
+# result[model]["LRT"] = chi2.sf(
+# result[model]["log(L)"] - result[basemodel]["log(L)"],
+# result[model]["k"] - result[basemodel]["k"],
+# )
+pd.DataFrame(result).T
 
 # %% [markdown]
-# \nocite{*}
-# \bibliographystyle{plain}
-# \bibliography{reference}
+# ## 第二题
+#
+# 对比五个模型的SBC和AIC数值，我们认为ARCH模型是最优的。因此，下文选择ARCH模型来估计投资组合的未来走势与市场风险。
+#
+# 我们使用VaR和ES来描述模拟分布的市场风险，分别探讨了置信水平为99%、95%和90%三种情况，最终得到使用GARCH模型进行蒙特卡洛模拟时，99%置信水平下的单周收益率VaR=-22.05%，对应的ES=-26.16%；95%置信水平下的单周收益率VaR=-15.17%，对应的ES=-19.54%；90%置信水平下的单周收益率VaR=-11.80%，对应的ES=-16.46%。
+
+# %%
+monte_carlo = {}
+np.random.RandomState(42)
+cases = 10000
+for case in range(cases):
+ if (case*10) % cases == 0:
+ print("simulating",case/cases*100, "%")
+ sim_model = arch_model(None, p=1, q=1)
+ simulation = sim_model.simulate(garch.params, 26)
+ monte_carlo[case] = np.prod(1 + simulation["data"])
+
+# %%
+monte_carlo = pd.Series(monte_carlo)
+var = pd.Series(monte_carlo).quantile([0.01, 0.05, 0.1])
+es = var.apply(lambda x: monte_carlo[monte_carlo < x].mean())
+pd.DataFrame({"es": es - 1, "var": var - 1})