Timing S&P500 Using a Large Set of Forecasting Variables

1from AlgorithmImports import *
2import data_tools
3from pandas.core.frame import DataFrame
4from dateutil.relativedelta import relativedelta
5import statsmodels.api as sm
6from statsmodels.tools.eval_measures import rmse
7# endregion
8class TimingSP500UsingaLargeSetofForecastingVariables(QCAlgorithm):
9def Initialize(self):
10self.SetStartDate(2000, 1, 1)
11self.SetCash(100000)
12self.period:int = 24
13self.quarter_period:int = 3 # due to Deliquencies quarterly data
14self.month_period:int = 21
15self.max_alocation:float = 1.5
16self.decay_factor:float = 0.99
17self.leverage:int = 3
18self.max_missing_days:int = 31
19self.independent_variables_num:int = 12
20self.market:Symbol = self.AddEquity("SPY", Resolution.Daily).Symbol
21self.bil:Symbol = self.AddEquity('BIL', Resolution.Daily).Symbol
22self.commodity_index:Symbol = self.AddEquity('GSG', Resolution.Daily).Symbol
23self.dollar_index:Symbol = self.AddEquity('UUP', Resolution.Daily).Symbol
24self.baltic_index:Symbol = self.AddData(data_tools.BalticDryIndex, 'BADI', Resolution.Daily).Symbol
25for symbol in [self.market, self.bil]:
26    self.Securities[symbol].SetLeverage(self.leverage)
27quandl_tickers:List[str] = ['RATEINF/CPI_USA']
28self.consumer_index_usa:Symbol = self.AddData(data_tools.QuandlValue, 'RATEINF/CPI_USA', Resolution.Daily).Symbol
29quarterly_custom_data:List[str] = ['DRALACBN_YOY']
30monthly_custom_data:List[str] = ['INDPRO_YOY', 'BAA_AAA', 'NEWORDER', 'HOUST', 'UNEMPLOYMENT_RATE']
31daily_custom_data:List[str] = ['T10Y3M']
32self.data:Dict[Symbol, float] = {}
33self.quarterly_custom_data:List[Symbol] = [self.AddData(data_tools.QuarterlyCustomData, x, Resolution.Daily).Symbol for x in quarterly_custom_data]
34self.monthly_custom_data:List[Symbol] = [self.AddData(data_tools.MonthlyCustomData, x, Resolution.Daily).Symbol for x in monthly_custom_data]
35self.daily_custom_data:List[Symbol] = [self.AddData(data_tools.DailyCustomData, x, Resolution.Daily).Symbol for x in daily_custom_data] + [self.baltic_index]
36self.qc_data:List[Symbol] = [self.consumer_index_usa]
37self.recent_month:int = -1
38def OnData(self, data: Slice):
39# monthly rebalance
40if self.recent_month == self.Time.month:
41    return
42self.recent_month = self.Time.month
43quarterly_custom_data_last_update_date:Dict[Symbol, datetime.date] = data_tools.QuarterlyCustomData._last_update_date
44monthly_custom_data_last_update_date:Dict[Symbol, datetime.date] = data_tools.MonthlyCustomData._last_update_date
45daily_custom_data_last_update_date:Dict[Symbol, datetime.date] = data_tools.DailyCustomData._last_update_date
46# check if data is still arriving
47if all([self.Securities[x].GetLastData() for x in self.quarterly_custom_data + self.monthly_custom_data + self.daily_custom_data + self.qc_data]) \
48    and any([self.Time.date() >= monthly_custom_data_last_update_date[x] for x in monthly_custom_data_last_update_date]) \
49    and any([self.Time.date() >= quarterly_custom_data_last_update_date[x] for x in quarterly_custom_data_last_update_date]) \
50    and any([self.Time.date() >= daily_custom_data_last_update_date[x] for x in daily_custom_data_last_update_date]) \
51    and any((self.Time.date() - self.Securities[x].GetLastData().Time.date()).days >= self.max_missing_days for x in self.qc_data):
52    self.Liquidate()
53    return Universe.Unchanged
54# call history on all variables
55history:DataFrame = self.History(self.monthly_custom_data + self.daily_custom_data + self.qc_data + self.quarterly_custom_data, start=self.Time.date() - relativedelta(months=self.period + self.quarter_period), end=self.Time.date())['value'].unstack(level=0)
56history_equities:DataFrame = self.History([self.market, self.commodity_index, self.dollar_index], start=self.Time.date() - relativedelta(months=self.period), end=self.Time.date())['close'].unstack(level=0)
57history = history.groupby(pd.Grouper(freq='MS')).last()
58history_equities = history_equities.groupby(pd.Grouper(freq='MS')).last()
59history['DRALACBN_YOY'] = history['DRALACBN_YOY'].fillna(method='ffill') # forward fill due quarterly data
60history = history[self.quarter_period:]
61history = pd.concat([history, history_equities], axis=1)
62history = history[1:]
63if history.iloc[0].isna().any():
64    return
65
66if len(history.columns) < self.independent_variables_num:
67    return
68history['normalized_market'] = history[self.market]
69history = history.fillna(method='ffill')
70history[['INDPRO_YOY', self.commodity_index, self.dollar_index, self.baltic_index, self.consumer_index_usa, 'normalized_market']] = history[['INDPRO_YOY', self.commodity_index, self.dollar_index, self.baltic_index, self.consumer_index_usa, 'normalized_market']].pct_change()
71history[['HOUST', 'UNEMPLOYMENT_RATE', 'DRALACBN_YOY']] = history[['HOUST', 'UNEMPLOYMENT_RATE', 'DRALACBN_YOY']].diff()
72# normalize variables with their 2 year standard deviation
73history_normalized = history.loc[:, history.columns != self.market].iloc[-1] / history.loc[:, history.columns != self.market][:-1].std()
74
75# save values to expanding list
76for symbol, value in history_normalized.items():
77    if symbol not in self.data:
78        self.data[symbol] = []
79    self.data[symbol].append(value)
80
81if self.market not in self.data:
82    self.data[self.market] = []
83self.data[self.market].append(history[self.market].pct_change().iloc[-1])
84    
85# run factor analysis
86if len(self.data[next(iter(self.data))]) >= self.period:
87    df:DataFrame = pd.DataFrame(self.data)
88    result_factors:Dict[Symbol, float] = self.factor_analysis(self.market, df, len(history.loc[:, history.columns != self.market].columns))
89    # run regression on best model and predict return of market 
90    if result_factors:
91        y:np.ndarray = df[self.market][1:].values
92        x:np.ndarray = df[list(result_factors.keys())][:-1].values
93        model = self.multiple_linear_regression(x, y)
94        predicted_y:np.ndarray = model.predict(sm.add_constant(df[list(result_factors.keys())][-1:].values, has_constant='add'))
95        rmse_value:float = rmse(y, predicted_y)
96        alocation:float = predicted_y * (1 / rmse_value) * 5
97        # trade execution
98        if predicted_y < 0:
99            self.Liquidate()
100            self.SetHoldings(self.bil, 1)
101            return
102        traded_alocation:Tuple[float, float] = (alocation if alocation < self.max_alocation else self.max_alocation, self.max_alocation - alocation if alocation < self.max_alocation else 0)
103        if self.market in data and data[self.market] and self.bil in data and data[self.bil]:
104            self.SetHoldings(self.market, traded_alocation[0])
105            self.SetHoldings(self.bil, traded_alocation[1])
106def factor_analysis(self, asset_id: str, factor_df: pd.DataFrame, n:int=0, ignore_negative_loading:bool = False) -> Dict:
107""" 
108Performs factor analysis.
109:param asset_id: Id of asset that is going to be fit.
110:param factor_df: Performance dataframe with all the factors and asset with asset_id.
111:param n: Max number of assets in the model.
112:return: Fitted model dictionary.
113"""
114order:List[str] = []
115y:np.ndarray = factor_df[asset_id].values 
116reduced_df:DataFrame = factor_df.drop([asset_id], axis=1)
117# reduced_df.drop(['Date'], axis=1, inplace=True)
118
119for i, _ in enumerate(reduced_df.items()):
120    # not negative beta models properties
121    aic_reduced:List[float] = []
122    col_names_reduced:List[str] = []
123        
124    for j, (j_asset_name, _) in enumerate(reduced_df.items()):
125        if j_asset_name in order:
126            continue
127        
128        if i > len(order) - 1:
129            order.append(j_asset_name)
130        else:
131            order[i] = j_asset_name
132        
133        order_reduced:List[str] = [x for x in order if x != '']
134        x:np.ndarray = reduced_df[order_reduced].values
135        model = self.multiple_linear_regression(x, y)
136        aic:float = model.aic
137        
138        if ignore_negative_loading:
139            # store only not negative beta models
140            if not any(x < 0. for x in betas_):
141                aic_reduced.append(aic)
142                col_names_reduced.append(j_asset_name)
143        else:
144            aic_reduced.append(aic)
145            col_names_reduced.append(j_asset_name)
146    
147    if len(aic_reduced) == 0:
148        order[i] = ''
149    else:
150        min_aic_index:int = min(range(len(aic_reduced)), key=aic_reduced.__getitem__)   # index of lowest aic
151        order[i] = col_names_reduced[min_aic_index]
152        
153    if i > 0:
154        order_reduced:List[str] = [x for x in order if x != '']
155        previous_order:List[str] = [x for x in order_reduced[:-1]]
156        
157        if len(order_reduced) != 0 and len(previous_order) != 0:
158            # build previous and actual regression models
159            x:np.ndarray = reduced_df[order_reduced].values
160            model = self.multiple_linear_regression(x, y)
161            aic:float = model.aic
162            betas:List[float] = model.params[1:]
163            # betas, _, lm_aic = self.multiple_linear_regression(x, y)
164            x = reduced_df[previous_order].values
165            model = self.multiple_linear_regression(x, y)
166            prev_aic:float = model.aic
167            prev_betas:List[float] = model.params[1:]
168            # prev_betas, _, previous_aic = self.multiple_linear_regression(y, x)
169            
170            if (aic >= prev_aic) or (n != 0 and len(order_reduced) > n):
171                result:Dict[str, float] = { order_id : beta for order_id, beta in zip(previous_order, betas.tolist()) }
172                return result
173
174order_reduced:List[str] = [x for x in order if x != '']
175x:np.ndarray = reduced_df[order_reduced].values
176model = self.multiple_linear_regression(x, y)
177betas:List[float] = model.params[1:]
178# betas, _, lm_aic = self.multiple_linear_regression(y, x)
179result:Dict[str, float] = { order_id : betas  for order_id, betas in zip(order_reduced, betas.tolist()) }
180return result
181def multiple_linear_regression(self, x:np.ndarray, y:np.ndarray):
182# x:np.ndarray = np.array(x).T
183x = sm.add_constant(x, prepend=True)
184w:List[float] = self.decay_factor ** np.arange(len(y), 0, -1)
185result = sm.WLS(endog=y, exog=x, weights=w).fit()
186return result