<jupyter_start><jupyter_text>0. ETF SelectionWe select the SPDR Gold Shares (GLD) ETF as the gold ETF. It is traded on Nasdaq, the currency is USD.Similarly, we choose the Amundi CAC 40 UCITS ETF-C (C40.PA) as the equity ETF. It will track the CAC 40 index of France. It is traded on Paris Euronext, the currency is EUR.The currency for Bitcoin is USD.Data source: https://finance.yahoo.com/  1. Data Importing<jupyter_code>import arch
import holidays
import pmdarima
import pandas as pd
import numpy as np
from pandas import Series, DataFrame
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from datetime import datetime
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.varmax import VARMAX
from statsmodels.tsa.stattools import adfuller, coint
from statsmodels.tsa.vector_ar.vecm import VECM
from sklearn import preprocessing
%matplotlib inline
gold_df = pd.read_csv("data/SPDR Gold Shares (GLD) Jan - Dec 2020.csv")
equity_df = pd.read_csv("data/Amundi CAC 40 UCITS ETF-C (C40.PA) Jan 2020 - Dec 2020.csv")
bitcoin_df = pd.read_csv('data/Bitcoin USD (BTC-USD) Jan 2020 - Dec 2020.csv')<jupyter_output><empty_output><jupyter_text>Convert the data into the datetime format and make it the index to query the dataframe easier.<jupyter_code>def convert_df(df):
    df["Date"] = pd.to_datetime(df["Date"], format="%Y-%m-%d")
    df.set_index("Date", inplace=True)   
    return df

gold_df = convert_df(gold_df)
equity_df = convert_df(equity_df)
bitcoin_df = convert_df(bitcoin_df)<jupyter_output><empty_output><jupyter_text>We use the common subset of days for all 3 time series to make them comparable and to run the test for cointegration in question 7.<jupyter_code>missing = set(equity_df.index) - set(gold_df.index)
def remove_missing_days(df, indices, missing):
    return df.loc[[index for index in indices if index not in missing]]
gold_df = remove_missing_days(gold_df, equity_df.index, missing)
equity_df = remove_missing_days(equity_df, equity_df.index, missing)
bitcoin_df = remove_missing_days(bitcoin_df, equity_df.index, missing)
assert equity_df.shape == gold_df.shape
assert equity_df.shape == bitcoin_df.shape<jupyter_output><empty_output><jupyter_text>We then fill some missing values in the equity series with its old values to avoid data leakage.<jupyter_code>equity_df = equity_df.ffill()<jupyter_output><empty_output><jupyter_text>Verify that the time range is correct.<jupyter_code>gold_df.head()
gold_df.tail()
equity_df.head()
equity_df.tail()
bitcoin_df.head()
bitcoin_df.tail()<jupyter_output><empty_output><jupyter_text>2. Data Processing We use adjusted close prices to calculate the daily returns. Adjusted close prices are the prices that already take into account stock split and dividends, which reflex more accurate the change of the prices.<jupyter_code>equity_df["Daily Return"] = equity_df["Adj Close"].pct_change(1)
equity_df.head()
gold_df["Daily Return"] = gold_df["Adj Close"].pct_change(1)
gold_df.head()
bitcoin_df["Daily Return"] = bitcoin_df["Adj Close"].pct_change(1)
bitcoin_df.head()<jupyter_output><empty_output><jupyter_text>5. Category 1 Models: Just use 1 variable.<jupyter_code># Helper
df_names = {0: "gold ETF", 1: "equity ETF", 2: "Bitcoin"}
dfs = [gold_df, equity_df, bitcoin_df]
def get_data(df, month_start, month_end, column=None):
    data = df[(df.index >= f"2020-{month_start:02d}-01") & (df.index < f"2020-{month_end:02d}-01")]
    if column: 
        data = data[column]
    return data
def summarize_data(df):
    data = get_data(df, 3, 12, "Adj Close")
    moving_avg = data.rolling(20, min_periods=1).mean()
    return get_data(moving_avg, 4, 12)

def get_data(df, month_start, month_end, column):
    return df[(df.index >= f"2020-{month_start:02d}-01") & (df.index < f"2020-{month_end:02d}-01")][column]
def fit_arima(data, exog= None):
    model = ARIMA(data, exog=exog, order=(2,0,2))
    model_fit = model.fit()
    return model_fit
def fit_garch(data, garch_type="GARCH"):
    if garch_type == "TARCH":
        garch = arch.arch_model(data, vol='TGARCH', p=1, o=1, q=1, power=1)
    else:
        garch = arch.arch_model(data, vol=garch_type, p=1, o=0, q=1)
    garch_fitted = garch.fit()
    print(garch_fitted.summary())
    if garch_type != "FIGARCH":
        omega = garch_fitted.params["omega"]
        alpha = garch_fitted.params["alpha[1]"]
        beta = garch_fitted.params["beta[1]"]
        print(f"Unconditional variance: {omega/(1 - alpha - beta)}")
    return garch_fitted
bitcoin_q2 = get_data(bitcoin_df, 4, 6, "Daily Return")
bitcoin_q3 = get_data(bitcoin_df, 7, 9, "Daily Return")
bitcoin_q4 = get_data(bitcoin_df, 10, 12, "Daily Return")
equity_q2 = get_data(equity_df, 4, 6, "Daily Return")
equity_q3 = get_data(equity_df, 7, 9, "Daily Return")
equity_q4 = get_data(equity_df, 10, 12, "Daily Return")
gold_q2 = get_data(gold_df, 4, 6, "Daily Return")
gold_q3 = get_data(gold_df, 7, 9, "Daily Return")
gold_q4 = get_data(gold_df, 10, 12, "Daily Return")
bitcoin_q2.shape
bitcoin_q3.shape
bitcoin_q4.shape
gold_q2.shape
gold_q3.shape<jupyter_output><empty_output><jupyter_text>Our first model with ARIMA is a simple one, we'll allocate our capital totally to bitcoin if we predict that the next day price will rise, otherwise we will short 100%.The total return for Q3 is 3.3%, not very impressive, but for Q4, it's 82%, which is much better. Return for 2 quarters is 87%<jupyter_code>model = fit_arima(bitcoin_q2)
return_pred = model.predict(start=bitcoin_q2.shape[0], end=bitcoin_q2.shape[0] + bitcoin_q3.shape[0] - 1).values
signal = np.where(return_pred > 0, 1, -1)
return_q3 = np.product(bitcoin_q3 * signal + 1) - 1
return_q3
model = fit_arima(bitcoin_q3)
return_pred = model.predict(start=bitcoin_q3.shape[0], end=bitcoin_q3.shape[0] + bitcoin_q4.shape[0] - 1).values
signal = np.where(return_pred >0, 1, -1)
return_q4 = np.product(bitcoin_q4 * signal + 1) - 1
return_q4
(1 + return_q3) * (1 + return_q4) - 1<jupyter_output><empty_output><jupyter_text>We prefer the ARMA model for the ease of model and a decent return  6. Category 2 Models: Just use 2 variables With VARMA model, we use the equity returns as to enhance the prediction for the bitcoin returns. for Q3, the return is 28% while for Q4, it's 63%.  Return for 2 quarters is 108%<jupyter_code>model = fit_arima(bitcoin_q2, equity_q2)
return_pred = model.predict(exog=equity_q3, start=bitcoin_q2.shape[0], end=bitcoin_q2.shape[0] + bitcoin_q3.shape[0] - 1).values
signal = np.where(return_pred > 0, 1, -1)
return_q3 = np.product(bitcoin_q3 * signal + 1) - 1
return_q3
model = fit_arima(bitcoin_q3, equity_q3)
return_pred = model.predict(exog=equity_q4, start=bitcoin_q3.shape[0], end=bitcoin_q3.shape[0] + bitcoin_q4.shape[0] - 1).values
signal = np.where(return_pred > 0, 1, -1)
return_q4 = np.product(bitcoin_q4 * signal + 1) - 1
return_q4
(1 + return_q3) * (1 + return_q4) - 1<jupyter_output><empty_output><jupyter_text>We prefer the VARMA model thanks to its execellent return  7 . Category 3 Models: Use all 3 variables  VARMA model With VARMA model, we use both the equity and gold returns as to enhance the prediction for the bitcoin returns. for Q3, the return is 50% while for Q4, it's 40%. Return for 2 quarters is 109%<jupyter_code>model = fit_arima(bitcoin_q2, pd.concat([equity_q2, gold_q2], axis=1))
return_pred = model.predict(exog=pd.concat([equity_q3, gold_q3], axis=1), start=bitcoin_q2.shape[0], end=bitcoin_q2.shape[0] + bitcoin_q3.shape[0] - 1).values
signal = np.where(return_pred > 0, 1, -1)
return_q3 = np.product(bitcoin_q3 * signal + 1) - 1
return_q3
model = fit_arima(bitcoin_q3, pd.concat([equity_q3, gold_q3], axis=1))
return_pred = model.predict(exog=pd.concat([equity_q4, gold_q4], axis=1), start=bitcoin_q3.shape[0], end=bitcoin_q3.shape[0] + bitcoin_q4.shape[0] - 1).values
signal = np.where(return_pred > 0, 1, -1)
return_q4 = np.product(bitcoin_q4 * signal + 1) - 1
return_q4
(1 + return_q3) * (1 + return_q4) - 1
# index = 2
# df = dfs[index]
# print(f"{model} model for {df_names[index]} from April to December")
# data = get_data(df, 4, 6, "Daily Return")
# data = data.dropna()
# garch_fitted = fit_garch(data, "GARCH")
# forecasts = garch_fitted.forecast(horizon=43)<jupyter_output><empty_output><jupyter_text>Compare models using cumulative return and volatility<jupyter_code>data1 = get_data(gold_df, 4, 6, "Daily Return").values
data2 = get_data(equity_df, 4, 6, "Daily Return").values
data3 = get_data(bitcoin_df, 4, 6, "Daily Return").values
coint(data1, np.array([data2, data3]).T)
coint(data2, np.array([data1, data3]).T)
coint(data3, np.array([data1, data2]).T)<jupyter_output><empty_output><jupyter_text>p-value < 0.01 for all tests, we conclude that there are cointegrating vectors.  Vector Error Correction Model We allocate our capital to 3 assets based on the predicted return, with a negative return means that we'll short the asset, and the sum of the 3 assets is always 1. Q3's return is 23%, Q4's return is 47%, total return is 81%<jupyter_code>data_q2 = np.array([data1, data2, data3]).T
model = VECM(data_q2, coint_rank=1)
vecm_res = model.fit()
res = vecm_res.predict(steps=43)
row_sums = res.sum(axis=1)
weights = res / row_sums[:, np.newaxis]
weights
data1 = get_data(gold_df, 7, 9, "Daily Return").values
data2 = get_data(equity_df, 7, 9, "Daily Return").values
data3 = get_data(bitcoin_df, 7, 9, "Daily Return").values
data_q3 = np.array([data1, data2, data3]).T
return_q3 = weights * data_q3
total_return_q3 = return_q3.sum(axis=1)
portfolio_return_q3 = np.product(total_return_q3 + 1) - 1
portfolio_return_q3
coint(data1, np.array([data2, data3]).T)
coint(data2, np.array([data1, data3]).T)
coint(data3, np.array([data1, data2]).T)<jupyter_output><empty_output><jupyter_text>Interestingly, we can't reject the null hypothesis for the test for the equity ETF and the combintion of gold and bitcoin ETF. We can reject the null hypothesis for the other two tests.<jupyter_code>model = VECM(data_q3, coint_rank=1)
vecm_res = model.fit()
res = vecm_res.predict(steps=42)
row_sums = res.sum(axis=1)
weights = res / row_sums[:, np.newaxis]
weights
data1 = get_data(gold_df, 10, 12, "Daily Return").values
data2 = get_data(equity_df, 10, 12, "Daily Return").values
data3 = get_data(bitcoin_df, 10, 12, "Daily Return").values
data_q4 = np.array([data1, data2, data3]).T
return_q4 = weights * data_q4
total_return_q4 = return_q4.sum(axis=1)
portfolio_return_q4 = np.product(total_return_q4 + 1) - 1
portfolio_return_q4
(1 + portfolio_return_q3) * (1 + portfolio_return_q4) - 1<jupyter_output><empty_output>