Net Emotional Volatility Index Effect

1from AlgorithmImports import *
2#endregion
3# https://quantpedia.com/strategies/net-emotional-volatility-index-effect/
4#
5# Historical S&P 500 index data are used to estimate the coefficients in the GARCH (1,1) model. Once the GARCH coefficients are obtained, for each day, the investor forecasts the volatility 
6# for the next 30 days, reiterating the forecast, which means that every forecast becomes the last observation for the next forecast. The average of these values is called the Average GARCH 
7# Forecast (AGF). The investor then calculates the Net Emotional Volatility Index (NEVI) as the difference between the VIX and AGF. The historical NEVI index is then sorted into deciles, and 
8# the current NEVI is used as a “Momentum” indicator – the investor opens long on a position when the indicator goes below the 2nd decile and short on positions when it goes above the 8th decile.
9# The position is held 30 days from the day of the last long or short signal.
10import numpy as np
11from arch import arch_model
12class NetEmotionalVolatilityIndexEffect(QCAlgorithm):
13def Initialize(self):
14self.SetStartDate(2008, 1, 1)
15self.SetCash(100000)
16
17self.period:int = 21
18self.holding_period:int = 30
19self.SetWarmUp(self.period + self.holding_period)
20    
21self.market:Symbol = self.AddEquity("SPY", Resolution.Daily).Symbol
22    
23self.vix:Symbol = self.AddData(QuandlVix, "CBOE/VIX", Resolution.Daily).Symbol # starts in 2004
24
25self.days_held:int = 0
26
27self.market_data:RollingWindow = RollingWindow[float](self.period)
28self.volatility:RollingWindow = RollingWindow[float](self.period)
29self.nevi_index:list[float] = []
30def OnData(self, data):
31# market data are present in the algorithm
32if self.market in data and data[self.market]:
33    market_price:float = data[self.market].Value
34    if market_price != 0:
35        # store market price data
36        self.market_data.Add(market_price)
37    
38    if self.market_data.IsReady:
39        # calcualte market daily returns
40        market_prices:np.ndarray = np.array([x for x in self.market_data][::-1])
41        market_returns:np.ndarray = (market_prices[:-1] / market_prices[1:] - 1) * 100
42        
43        # forecast market volatility
44        exp_market_volatility:float = self.garch_vol_predict_volatility(market_returns, 1, 1, 30)
45        self.volatility.Add(exp_market_volatility)
46        
47        if self.volatility.IsReady:
48            agf:float = np.mean([x for x in self.volatility])
49            
50            # vix data are present in the algorithm
51            if self.vix in data and data[self.vix]:
52                vix:float = data[self.vix].Value
53                self.Plot("Expected Market Volatility", "Expected Market Volatility", exp_market_volatility)
54                self.Plot("VIX", "VIX", vix)
55                if vix != 0:
56                    nevi:float = vix - agf
57                    self.nevi_index.append(nevi)
58                    
59                    # wait until at least 30 data points is present for NEVI
60                    if len(self.nevi_index) >= self.period and not self.Portfolio[self.market].Invested and not self.IsWarmingUp:
61                        low_decile:float = np.percentile(self.nevi_index[:-1], 20)
62                        high_decile:float = np.percentile(self.nevi_index[:-1], 80)
63                        
64                        recent_nevi:float = self.nevi_index[-1]
65                        if recent_nevi > high_decile:
66                            self.SetHoldings(self.market, -1)
67                            self.days_held = 0
68                        elif recent_nevi < low_decile:
69                            self.SetHoldings(self.market, 1)
70                            self.days_held = 0
71
72if self.Portfolio[self.market].Invested:
73    self.days_held += 1
74    if self.days_held == self.holding_period:
75        self.Liquidate(self.market)
76        self.days_held = 0
77def garch_vol_predict_volatility(self, returns, p, q, h) -> float:
78am = arch_model(returns, mean='Zero', vol='GARCH', p=p, q=q)
79res = am.fit(disp='off')
80forecasts = res.forecast(horizon=h)
81return np.sqrt(forecasts.variance.values[-1,:][0] * 252)
82
83class QuandlVix(PythonQuandl):
84def __init__(self):
85self.ValueColumnName = "close"