Forecasting the Nasdaq-100's volatility with HAR-RV

1"""A HAR-RV volatility forecast as a filter for the 15-minute opening range breakout, on Micro E-mini Nasdaq-100 (MNQ).2 3The forecast is Corsi's (2009) HAR model of realized volatility, in logs, fitted day by day on the past4only (an expanding window), so every forecast is out of sample when it is made.5 6Rules:7  1. Realized volatility (RV) of a session is the square root of the sum of its squared 5-minute log8     returns, regular session only (the overnight gap left out).9  2. After each session, the model log RV(t+1) = b0 + b1 log RV(t) + b2 mean(log RV, last 5 days) + b310     mean(log RV, last 22 days) is refitted by least squares on every past day, and forecasts the next.11  3. Trade the next day only if its forecast is at least the `min_quantile` (0.5, the median) of the12     last 250 forecasts, and at most their `max_quantile` (1), and after 250 days of forecasts.13  4. The trade is the 15-minute opening range breakout: the first 5-minute close outside the 09:30-09:4514     ET range, until 12:00 ET, stop at the other side, target `target_r` (2) times the risk, flat at15     15:55 ET.16  5. Size: about `risk_usd` dollars of risk a trade (micro contracts, rounded down, at least 1 and at17     most `max_contracts`).18"""19import math20 21import qy22 23 24def ready(*values):25    """True when every value is a number (indicators have no value until they have enough data)."""26    return all(v is not None and v == v for v in values)27 28 29def solve(a, b):30    """Gaussian elimination with partial pivoting for a small system a x = b."""31    n = len(b)32    m = [row[:] + [b[i]] for i, row in enumerate(a)]33    for c in range(n):34        p = max(range(c, n), key=lambda r: abs(m[r][c]))35        if abs(m[p][c]) < 1e-12:36            return None37        m[c], m[p] = m[p], m[c]38        for r in range(n):39            if r != c:40                f = m[r][c] / m[c][c]41                m[r] = [x - f * y for x, y in zip(m[r], m[c])]42    return [m[i][n] / m[i][i] for i in range(n)]43 44 45class Bot(qy.Strategy):46    """Trade the opening range breakout only on days the HAR model expects to be volatile."""47 48    symbol = "MNQ"49    timeframe = "5m"50    session = "rth"51    params = {"min_quantile": 0.5, "max_quantile": 1.0, "target_r": 2.0, "risk_usd": 250, "max_contracts": 50}52    max_position = 5053    entry_window = ("09:45", "12:00")54    flatten_time = "15:55"55    overnight = False56 57    def setup(self):58        self.or_high, self.or_low = self.opening_range(minutes=15)59        self.day = None60        self.last = None61        self.ss = 0.062        self.logrv = []63        self.xtx = [[0.0] * 4 for _ in range(4)]64        self.xty = [0.0] * 465        self.n = 066        self.forecasts = []67        self.trade_today = False68        self.done_day = None69 70    def _size(self, risk_points):71        per_contract = risk_points * self.point_value72        n = int(self.p.risk_usd // per_contract) if per_contract > 0 else 173        return max(1, min(self.p.max_contracts, n))74 75    def _features(self, k):76        lr = self.logrv77        return [1.0, lr[k], sum(lr[k - 4:k + 1]) / 5, sum(lr[k - 21:k + 1]) / 22]78 79    def _end_of_day(self):80        if self.ss <= 0:81            return82        self.logrv.append(math.log(math.sqrt(self.ss) * 100))83        k = len(self.logrv) - 184        if k >= 22:  # yesterday's features, today's value: one more pair for the fit85            x = self._features(k - 1)86            for i in range(4):87                self.xty[i] += x[i] * self.logrv[k]88                for j in range(4):89                    self.xtx[i][j] += x[i] * x[j]90            self.n += 191        if k >= 21 and self.n >= 60:92            beta = solve(self.xtx, self.xty)93            if beta:94                self.forecasts.append(sum(b * v for b, v in zip(beta, self._features(k))))95                self.forecasts = self.forecasts[-251:]96 97    def _decide(self):98        past = self.forecasts[:-1]99        if len(past) < 250:100            return False101        f = self.forecasts[-1]102        rank = sum(1 for v in past if v <= f) / len(past)103        return self.p.min_quantile <= rank <= self.p.max_quantile104 105    def on_bar(self, bar):106        if bar.date != self.day:107            if self.day is not None:108                self._end_of_day()109            self.day, self.ss, self.last = bar.date, 0.0, bar.close110            self.trade_today = self._decide()111        else:112            r = math.log(bar.close / self.last)113            self.ss += r * r114            self.last = bar.close115        high, low = self.or_high[0], self.or_low[0]116        if not self.trade_today or not ready(high, low) or self.done_day == bar.date or not self.position.is_flat:117            return118        close = bar.close119        if low <= close <= high:120            return121        self.done_day = bar.date122        if bar.dt.hour * 60 + bar.dt.minute + 5 > 720:123            return124        f = (high - low) / close * 100125        if close > high:126            self.buy(self._size(close - low), stop_price=low, target_r=self.p.target_r, tag=f"HAR ORB long | f={f:.4f}")127        else:128            self.sell(self._size(high - close), stop_price=high, target_r=self.p.target_r, tag=f"HAR ORB short | f={f:.4f}")