前言

美赛已经结束4天了，一直忙于教资考试的准备，今天我终于抽空写了这篇C题思路复盘的博客。

题目大致要求

题目叫'Trading Strategies（交易策略）'，一共给了两个文件，分别是比特币和黄金价格随时间变化的CSV文件。大致要求可以分为以下四个步骤：

基于截止至当日的价格情况建立模型，预测2021年9月10日原来的本金1000美元会变成多少钱？（Develop a model that gives the best daily trading strategy based only on price data up to that day. How much is the initial $1000 investment worth on 9/10/2021 using your model and strategy?）
证明你的模型提供了最佳策略（Present evidence that your model provides the best strategy.）
就是让你对手续费那两个参数做敏感性分析（Determine how sensitive the strategy is to transaction costs. How do transaction costs affect the strategy and results?）
美赛的老一套，让你写一个什么备忘录，向投资者陈述你的战略、模型和结果（Communicate your strategy, model, and results to the trader in a memorandum of at most two pages.）

具体题目内容见美赛网站：2022 MCM Problem C (immchallenge.org)http://www.immchallenge.org/mcm/2022_MCM_Problem_C.pdf

分析题目?

对于问题1，思路无非是分为两个步骤，第一个步骤是预测出后几日的价格进而估计得到后几日的收益率，然后通过动态规划模型进行最优化的求解。具体思路如下：

首先先看数据集，你会发现数据集只有两列，时间一列和价格一列。这意味着知网或者其他数据库中查到的很多论文写的预测算法你都用不了了。为什么呢？因为人家写论文时候用的数据集都有其他的特征，而不是这里只有时间一列特征。所以那些什么SVM、贝叶斯网、向量自回归等要求多维特征的算法都不适用。查来查去，有推荐使用LSTM、神经网络模型还有xgboost算法，但我都是不是很熟悉，最终还是用了时间序列ARIMA模型来做的预测。为什么用ARIMA模型来做的预测，我给出以下几点原因：

ARIMA模型属于单时间序列模型，符合题目数据集提供的数据以及题目给出的要求。
ARIMA模型是一个很常见的时间序列模型，较为成熟也是我所熟悉的，我能拿它写不少东西。
我使用Python编程，而Python的statsmodels库中提供了ARIMA现成的API，不用另外写代码，而且调用API绘制图像的操作十分方便。

对于问题二呢，网上的思路都是说，改变什么参数，看看收益率会不会提高，如果不会的话就说明是最优的，但也没说清楚具体怎么操作。而我的想法不完全是这样的，我看了他们的思路突然想到了"梯度下降法"，那你说如果我用模型的两个参数（具体来讲一个是ARIMA模型的AR参数另一个是MA参数，作为梯度下降法改变模型的两个搜索方向，然后进行算法计算，理论上是完全可行的，于是乎我就这么干了。

对于问题三和问题四就不再多说了。但我也罗嗦一下，问敏感性分析，我除了做了问题三的手续费的敏感性分析外，我还针对ARIMA模型做了敏感度分析。问题四就总分总来写，总一段，分三段，也就是战略、模型和结果各一段写，最后写个总结，要注意的是英文备忘录是有格式的，得按照格式来，要有点仪式感，不然你就可能会被扣分。

还有就是别忘了一开始写数据清洗那一段，我用了基于时间序列的牛顿插值法进行插补，用箱型图法判断离群点。

接下来我贴出我的代码。

代码

import pandas as pd
from pylab import *
import statsmodels.api as sm
from datetime import datetime
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.stattools import kpss
from statsmodels.stats.diagnostic import acorr_ljungbox


def adf_test(timeseries):
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value', '#Lags Used', 'Number of Observations Used'])
    for key, value in dftest[4].items():
       dfoutput['Critical Value (%s)'%key] = value
    print(dfoutput)


def kpss_test(timeseries):
    print('Results of KPSS Test:')
    kpsstest = kpss(timeseries, regression='c')
    kpss_output = pd.Series(kpsstest[0:3], index=['Test Statistic','p-value','Lags Used'])
    for key, value in kpsstest[3].items():
        kpss_output['Critical Value (%s)'%key] = value
    print(kpss_output)


dataGOLD = pd.read_csv("C:\\ProblemC\\LBMA-GOLD.csv")
dataMKPRU = pd.read_csv("C:\\ProblemC\\BCHAIN-MKPRU.csv")
dataGOLD['Date'] = pd.to_datetime(dataGOLD.Date)
dataMKPRU['Date'] = pd.to_datetime(dataMKPRU.Date)
dataGOLD = dataGOLD[dataGOLD['Date'] >= datetime.strptime('2016-9-11', "%Y-%m-%d")]
dataMKPRU = dataMKPRU[dataMKPRU['Date'] >= datetime.strptime('2016-9-11', "%Y-%m-%d")]
dataGOLD = dataGOLD[dataGOLD['Date'] <= datetime.strptime('2021-9-10', "%Y-%m-%d")]
dataMKPRU = dataMKPRU[dataMKPRU['Date'] <= datetime.strptime('2021-9-10', "%Y-%m-%d")]
dataGOLD = dataGOLD.sort_values('Date').reset_index(drop=True)
dataMKPRU = dataMKPRU.sort_values('Date').reset_index(drop=True)

adf_test(dataGOLD['Value'])
kpss_test(dataGOLD['Value'])

diffGOLD = dataGOLD['Value'].diff(1).dropna()
diffMKPRU = dataMKPRU['Value'].diff(1).dropna()

fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(diffMKPRU, ax=ax1)
ax1.xaxis.set_ticks_position('top')
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(diffMKPRU, ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
fig.tight_layout()
plt.show()

fig = plt.figure(figsize=(12, 8))
ax11 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(diffGOLD, ax=ax11)
ax11.xaxis.set_ticks_position('top')
ax22 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(diffGOLD, ax=ax22)
ax22.xaxis.set_ticks_position('bottom')
fig.tight_layout()
plt.show()

p_value = acorr_ljungbox(dataGOLD['Value'])
print(p_value)

results = sm.tsa.arma_order_select_ic(dataMKPRU['Value'], ic=['aic', 'bic'], trend='nc', max_ar=5, max_ma=5)
print('AIC', results.aic_min_order)
print('BIC', results.bic_min_order)

n_sample = dataGOLD['Value'].shape[0]
n_train = int(n_sample * 0.95) + 1
n_forecast = n_sample - n_train

ts_train = dataGOLD.iloc[:n_train]['Value']
ts_test = dataGOLD.iloc[n_train:]['Value']

arima = sm.tsa.SARIMAX(ts_train, order=(2, 1, 0))
model_results = arima.fit()
model_results.plot_diagnostics(figsize=(16, 12))
plt.show()
print(model_results.summary())

plt.title("Bitcoin price image after first-order difference")
plt.plot(dataMKPRU['Date'], dataMKPRU['Value'].diff(1))
plt.show()

plt.title("Gold price image after first-order difference")
plt.plot(dataGOLD['Date'], dataGOLD['Value'].diff(1))
plt.show()