import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"
import numpy as np, pandas as pd, matplotlib.pyplot as plt, time, warnings
warnings.filterwarnings("ignore")
import torch, torch.nn as nn; torch.set_num_threads(2)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("spx_rv_ret.csv"); dates=pd.to_datetime(d["date"]); rv=d["rv"].values; ret=d["ret"].values; n=len(rv)
lvol=np.log(np.sqrt(rv)*100.0)                        # log realized vol (%) -- the target
acf1=lambda x: np.corrcoef(x[1:],x[:-1])[0,1]
print(f"S&P realized volatility: {n} days, {d['date'].iloc[0]} to {d['date'].iloc[-1]}")
print(f"ACF(1)  log-RV = {acf1(lvol):.3f}  (persistent -> predictable)   returns = {acf1(ret):.3f}  (~0 -> unpredictable)")
fig,ax=plt.subplots(1,2,figsize=(14,4))
ax[0].plot(dates,np.sqrt(rv)*100*np.sqrt(252),color=BLUE,lw=.7); ax[0].set_title("S&P realized volatility (annualized %) — clustering, 2008 spike"); ax[0].set_ylabel("annualized vol %")
lags=range(1,41)
ax[1].bar([l-0.2 for l in lags],[np.corrcoef(lvol[l:],lvol[:-l])[0,1] for l in lags],width=.4,color=BLUE,label="log-RV")
ax[1].bar([l+0.2 for l in lags],[np.corrcoef(ret[l:],ret[:-l])[0,1] for l in lags],width=.4,color=GREY,label="returns")
ax[1].axhline(0,color="k",lw=.5); ax[1].set_xlabel("lag (days)"); ax[1].set_ylabel("autocorrelation"); ax[1].set_title("Volatility is autocorrelated; returns are not"); ax[1].legend()
plt.tight_layout(); plt.show()
print("The right panel is the whole motivation: volatility has long, slowly-decaying memory (forecastable); returns have")
print("none (a random walk). So we forecast volatility -- and a sequence model should be able to exploit that memory.")

S&P realized volatility: 3459 days, 2000-01-03 to 2013-11-12
ACF(1)  log-RV = 0.784  (persistent -> predictable)   returns = -0.087  (~0 -> unpredictable)

The right panel is the whole motivation: volatility has long, slowly-decaying memory (forecastable); returns have
none (a random walk). So we forecast volatility -- and a sequence model should be able to exploit that memory.

def rnn_forward(x_seq, Wx, Wh, bx, bh):
    H=Wh.shape[0]; h=np.zeros(H); out=[]
    for t in range(len(x_seq)):
        h=np.tanh(Wx@x_seq[t] + bx + Wh@h + bh)      # same weights every step (weight sharing across time)
        out.append(h)
    return np.array(out)
torch.manual_seed(0)
rnn=nn.RNN(input_size=1,hidden_size=8,batch_first=True)
xseq=np.random.default_rng(0).standard_normal((30,1)).astype(np.float32)
with torch.no_grad(): ref,_=rnn(torch.tensor(xseq).unsqueeze(0)); ref=ref[0].numpy()
p=dict(rnn.named_parameters())
mine=rnn_forward(xseq, p["weight_ih_l0"].detach().numpy(), p["weight_hh_l0"].detach().numpy(),
                 p["bias_ih_l0"].detach().numpy(), p["bias_hh_l0"].detach().numpy())
print("from-scratch RNN vs torch nn.RNN: max|diff| of hidden states =", float(np.abs(mine-ref).max()))
print("Identical -- nn.RNN is exactly this tanh recurrence with weights shared across all time steps.")
print("\nThe catch: unrolled over T steps, the gradient is multiplied by W_h T times when it flows back. If its eigenvalues")
print("are <1 the signal vanishes (long-range memory is lost); >1 and it explodes. That vanishing-gradient problem is what")
print("the LSTM was invented to solve.")

from-scratch RNN vs torch nn.RNN: max|diff| of hidden states = 8.327705092181503e-08
Identical -- nn.RNN is exactly this tanh recurrence with weights shared across all time steps.

The catch: unrolled over T steps, the gradient is multiplied by W_h T times when it flows back. If its eigenvalues
are <1 the signal vanishes (long-range memory is lost); >1 and it explodes. That vanishing-gradient problem is what
the LSTM was invented to solve.

L=22
# The window ends at t-1 and the target is t: a genuine ONE-step-ahead forecast, on the same information
# set HAR, AR(1) and the random walk are given. Targeting lvol[t+1] from a window ending at t-1 would skip
# day t entirely -- the single most informative predictor, correlated 0.78 with the target -- and quietly
# turn this into a two-step-ahead problem, handicapping the network against its own benchmarks.
Xs=np.array([lvol[t-L:t] for t in range(L,n)]); ys=np.array([lvol[t] for t in range(L,n)])
sp=int(0.8*len(ys))                                   # train 2000-~2010, test ~2011-2013
mu,sd=Xs[:sp].mean(),Xs[:sp].std()
Xt=torch.tensor(((Xs-mu)/sd).reshape(-1,L,1),dtype=torch.float32); yt=torch.tensor((ys-mu)/sd,dtype=torch.float32).view(-1,1)
class LSTMForecaster(nn.Module):
    def __init__(self,hidden=32):
        super().__init__(); self.lstm=nn.LSTM(1,hidden,batch_first=True); self.head=nn.Linear(hidden,1)
    def forward(self,x): out,_=self.lstm(x); return self.head(out[:,-1])
torch.manual_seed(0); lstm=LSTMForecaster(32)
opt=torch.optim.Adam(lstm.parameters(),lr=5e-3); lf=nn.MSELoss(); t=time.time()
for ep in range(40):
    perm=torch.randperm(sp)
    for b0 in range(0,sp,64):
        bi=perm[b0:b0+64]; opt.zero_grad(); lf(lstm(Xt[bi]),yt[bi]).backward(); opt.step()
with torch.no_grad(): lstm_pred=lstm(Xt[sp:]).numpy().ravel()*sd+mu
print(f"LSTM trained in {time.time()-t:.0f}s;  {sp} train windows / {len(ys)-sp} test windows (22-day input)")
print(f"LSTM params: {sum(p.numel() for p in lstm.parameters())}")

LSTM trained in 3s;  2749 train windows / 688 test windows (22-day input)
LSTM params: 4513

from numpy.linalg import lstsq
from arch import arch_model
def rmse(a,b): return np.sqrt(np.mean((a-b)**2))
day=lvol; wk=pd.Series(lvol).rolling(5).mean().values; mo=pd.Series(lvol).rolling(22).mean().values
Xh=np.column_stack([day,wk,mo]); yh=np.r_[lvol[1:],np.nan]
ok=~np.isnan(Xh).any(1)&~np.isnan(yh); Xh,yhv=Xh[ok],yh[ok]; s=int(0.8*len(yhv))
res={}
res["random walk"]=rmse(Xh[s:,0],yhv[s:])
ba=lstsq(np.column_stack([np.ones(s),Xh[:s,0]]),yhv[:s],rcond=None)[0]; res["AR(1)"]=rmse(ba[0]+ba[1]*Xh[s:,0],yhv[s:])
bh=lstsq(np.column_stack([np.ones(s),Xh[:s]]),yhv[:s],rcond=None)[0]; res["HAR-RV"]=rmse(np.column_stack([np.ones(len(Xh)-s),Xh[s:]])@bh,yhv[s:])
split=int(0.8*n); gr=arch_model(ret,mean="Constant",vol="GARCH",p=1,q=1,dist="t").fit(last_obs=split,disp="off")
fc=gr.forecast(horizon=1,start=split,reindex=False); gv=np.log(np.sqrt(fc.variance.values[:,0])); tg=lvol[split+1:split+1+len(gv)]; mm=min(len(gv),len(tg))
res["GARCH(1,1)-t"]=rmse(gv[:mm],tg[:mm])
res["LSTM"]=rmse(lstm_pred,ys[sp:])
tab=pd.Series(res).sort_values()
print("OOS RMSE on log realized volatility (lower = better):"); print(tab.round(4).to_string())
fig,ax=plt.subplots(figsize=(8,4.2)); cols=[BLUE if k=="LSTM" else (GREEN if k=="HAR-RV" else GREY) for k in tab.index]
ax.barh(tab.index,tab.values,color=cols); ax.invert_yaxis(); ax.set_xlabel("OOS RMSE (log realized vol)"); ax.set_title("Volatility forecast horse race")
for i,v in enumerate(tab.values): ax.text(v+0.002,i,f"{v:.3f}",va="center",fontsize=9)
plt.tight_layout(); plt.show()
har_pred=np.column_stack([np.ones(len(Xh)-s),Xh[s:]])@bh
_act=ys[sp:]
assert len(har_pred)==len(lstm_pred)==len(_act) and np.allclose(_act, yhv[s:]), "test sets must align to compare"
def _nw(x,lag=5):
    x=x-x.mean(); v=np.mean(x*x)
    for k in range(1,lag+1): v+=2*(1-k/(lag+1))*np.mean(x[k:]*x[:-k])
    return np.sqrt(v/len(x))
_dl=(lstm_pred-_act)**2-(har_pred-_act)**2; _dm=_dl.mean()/_nw(_dl)
_bias=float(np.mean(gv[:mm]-tg[:mm])); _gbc=rmse(gv[:mm]-_bias,tg[:mm])
print(f"\nHAR-RV {res['HAR-RV']:.4f} and the LSTM {res['LSTM']:.4f} are separated by {abs(res['HAR-RV']-res['LSTM']):.4f} RMSE -- which is not a")
print(f"result until it is tested. Both forecast the SAME {len(_act)} test days from the same information set, so a")
print(f"Diebold-Mariano test on squared-error loss applies directly: DM = {_dm:+.2f} against a 1.96 critical value.")
print("They are statistically indistinguishable. Neither model wins; the interesting fact is that they tie.")
print(f"\nThat tie is the finding worth carrying: HAR-RV spends {len(bh)} coefficients where the LSTM spends")
print(f"{sum(p.numel() for p in lstm.parameters()):,} -- about {sum(p.numel() for p in lstm.parameters())//len(bh):,} times as many parameters -- to reach the same accuracy. Three")
print("hand-designed features (yesterday, last week, last month) capture essentially everything a recurrent network")
print("can extract from this series on its own history.")
# A bias estimated on the test set is an ORACLE quantity -- it uses the answers. To claim the offset is a
# real, removable feature rather than a curve-fit, re-estimate it on the TRAINING period only and carry it
# forward; that is a correction someone could actually have applied in advance.
_gin=gr.conditional_volatility[:split]; _okv=~np.isnan(_gin)
_btr=float(np.mean(np.log(_gin[_okv])-lvol[:split][_okv])); _gtr=rmse(gv[:mm]-_btr,tg[:mm])
_overnight=1-(rv*1e4).mean()/ret.var(); _lognight=float(np.log(np.sqrt(ret.var()/(rv*1e4).mean())))
print(f"\nBoth clear AR(1) ({res['AR(1)']:.4f}) and the random walk ({res['random walk']:.4f}). GARCH ({res['GARCH(1,1)-t']:.4f}) trails -- but not")
print("for the reason usually given, and the arithmetic is worth doing rather than asserting.")
print(f"\nIts forecasts sit {_bias:+.3f} above realized log-vol on average, and that single constant accounts for")
print(f"{100*_bias**2/np.mean((gv[:mm]-tg[:mm])**2):.0f}% of its squared error. Two questions follow. Is the offset REAL or just fitted? Re-estimating it on")
print(f"the training period alone gives {_btr:+.3f} -- close to the {_bias:+.3f} the test set implies -- and applying that")
print(f"forward-looking-free correction gives RMSE {_gtr:.4f}, against {_gbc:.4f} for the oracle version. The level shift is")
print("stable, so this is a genuine property of the comparison and not an artefact of peeking at the answers.")
print(f"\nAnd WHERE does it come from? Partly the overnight gap: realized volatility is built from intraday returns")
print(f"only, while GARCH models close-to-close ones. Overnight moves are {100*_overnight:.0f}% of total daily variance here,")
print(f"which on its own implies a log offset of {_lognight:+.3f} -- about half of the {_bias:+.3f} observed. The remainder is the")
print(f"plain unconditional level gap between the two series: sd(returns) = {ret.std():.2f} against mean realized vol")
print(f"{(np.sqrt(rv)*100).mean():.2f}, a ratio of {ret.std()/(np.sqrt(rv)*100).mean():.2f}. GARCH is not forecasting badly; it is forecasting a different quantity.")
print(f"\nCorrected, GARCH ({_gtr:.4f}) beats the random walk ({res['random walk']:.4f}) and still trails HAR ({res['HAR-RV']:.4f}), with a")
print(f"correlation of {np.corrcoef(gv[:mm],tg[:mm])[0,1]:.2f} against realized vol. That ordering -- GARCH ahead of the random walk -- is the one")
print("the Realized Volatility notebooks report on their own scale, so the two projects agree once the level is")
print("reconciled. The raw ranking here was the odd one out, and the offset is why.")

OOS RMSE on log realized volatility (lower = better):
LSTM            0.3465
HAR-RV          0.3473
AR(1)           0.3916
random walk     0.4206
GARCH(1,1)-t    0.4568

HAR-RV 0.3473 and the LSTM 0.3465 are separated by 0.0009 RMSE -- which is not a
result until it is tested. Both forecast the SAME 688 test days from the same information set, so a
Diebold-Mariano test on squared-error loss applies directly: DM = -0.35 against a 1.96 critical value.
They are statistically indistinguishable. Neither model wins; the interesting fact is that they tie.

That tie is the finding worth carrying: HAR-RV spends 4 coefficients where the LSTM spends
4,513 -- about 1,128 times as many parameters -- to reach the same accuracy. Three
hand-designed features (yesterday, last week, last month) capture essentially everything a recurrent network
can extract from this series on its own history.

Both clear AR(1) (0.3916) and the random walk (0.4206). GARCH (0.4568) trails -- but not
for the reason usually given, and the arithmetic is worth doing rather than asserting.

Its forecasts sit +0.259 above realized log-vol on average, and that single constant accounts for
32% of its squared error. Two questions follow. Is the offset REAL or just fitted? Re-estimating it on
the training period alone gives +0.221 -- close to the +0.259 the test set implies -- and applying that
forward-looking-free correction gives RMSE 0.3783, against 0.3764 for the oracle version. The level shift is
stable, so this is a genuine property of the comparison and not an artefact of peeking at the answers.

And WHERE does it come from? Partly the overnight gap: realized volatility is built from intraday returns
only, while GARCH models close-to-close ones. Overnight moves are 23% of total daily variance here,
which on its own implies a log offset of +0.128 -- about half of the +0.259 observed. The remainder is the
plain unconditional level gap between the two series: sd(returns) = 1.32 against mean realized vol
0.97, a ratio of 1.36. GARCH is not forecasting badly; it is forecasting a different quantity.

Corrected, GARCH (0.3783) beats the random walk (0.4206) and still trails HAR (0.3473), with a
correlation of 0.66 against realized vol. That ordering -- GARCH ahead of the random walk -- is the one
the Realized Volatility notebooks report on their own scale, so the two projects agree once the level is
reconciled. The raw ranking here was the odd one out, and the offset is why.

# ys[i] is lvol at index L+i, so the test targets are dates[L:][sp:]. Both models now forecast the
# SAME days, so one date vector serves both -- which is also what makes the paired test above valid.
test_dates=dates.values[L:][sp:]
assert len(test_dates)==len(lstm_pred)==len(har_pred)
fig,ax=plt.subplots(figsize=(12,4.2))
ax.plot(test_dates, np.exp(ys[sp:]), color="black", lw=1.1, label="realized volatility (actual)")
ax.plot(test_dates, np.exp(lstm_pred), color=BLUE, lw=1.1, alpha=.9, label=f"LSTM (RMSE {res['LSTM']:.3f})")
ax.plot(test_dates, np.exp(har_pred), color=GREEN, lw=1.1, alpha=.8, label=f"HAR-RV (RMSE {res['HAR-RV']:.3f})")
ax.set_ylabel("realized volatility (%, daily)"); ax.set_title("Out-of-sample volatility forecasts vs actual (2011-2013)"); ax.legend()
plt.tight_layout(); plt.show()
print("Both models track the volatility dynamics closely -- rising into turbulent stretches, falling in calm ones. HAR is")
print("slightly sharper; the LSTM occasionally over-smooths. Neither is fooled by the persistence the random walk relies on.")

Both models track the volatility dynamics closely -- rising into turbulent stretches, falling in calm ones. HAR is
slightly sharper; the LSTM occasionally over-smooths. Neither is fooled by the persistence the random walk relies on.

fig,ax=plt.subplots(1,2,figsize=(12,4.6))
for a,(nm,pred,act,c) in zip(ax,[("LSTM",lstm_pred,ys[sp:],BLUE),("HAR-RV",har_pred,yhv[s:],GREEN)]):
    assert len(pred)==len(act)
    a.scatter(pred,act,s=8,alpha=.25,color=GREY); lim=[min(pred.min(),act.min()),max(pred.max(),act.max())]
    qd=np.quantile(pred,np.linspace(0,1,11)); bd=np.clip(np.digitize(pred,qd[1:-1]),0,9)
    pmk=[pred[bd==k].mean() for k in range(10)]; amk=[act[bd==k].mean() for k in range(10)]
    a.plot(lim,lim,"k--",lw=1,label="perfect"); a.plot(pmk,amk,"o-",color=c,lw=2,label="decile means")
    a.set_xlabel(f"{nm} predicted log-vol"); a.set_ylabel("actual log-vol")
    a.set_title(f"{nm}: proportions vs predictions (RMSE {rmse(pred,act):.3f})"); a.legend(fontsize=8)
plt.tight_layout(); plt.show()
# "no drift" is a claim about the decile means, so measure it rather than eyeball the cloud
for nm,pred,act in [("LSTM",lstm_pred,ys[sp:]),("HAR-RV",har_pred,yhv[s:])]:
    qd=np.quantile(pred,np.linspace(0,1,11)); bd=np.clip(np.digitize(pred,qd[1:-1]),0,9)
    g=[act[bd==k].mean()-pred[bd==k].mean() for k in range(10)]
    sl=np.polyfit(pred,act,1)
    print(f"{nm:7s} decile gaps run {min(g):+.3f} to {max(g):+.3f} (mean {np.mean(g):+.3f}); "
          f"regression of actual on predicted has slope {sl[0]:.3f}")
print("\nA scatter cloud cannot show drift; the decile means can. Both forecasters track the diagonal closely, with")
print("gaps small relative to the 0.35 RMSE and slopes near the 1.0 that an unbiased forecast requires -- so the")
print("'well-calibrated, no systematic over- or under-prediction' reading holds up when it is actually measured.")
print("That matters more than RMSE for risk use: a volatility forecast biased high or low mis-sizes every position")
print("built on it, and the bias would be invisible in the error metric alone.")

LSTM    decile gaps run -0.090 to +0.074 (mean -0.008); regression of actual on predicted has slope 0.954
HAR-RV  decile gaps run -0.100 to +0.060 (mean -0.007); regression of actual on predicted has slope 0.929

A scatter cloud cannot show drift; the decile means can. Both forecasters track the diagonal closely, with
gaps small relative to the 0.35 RMSE and slopes near the 1.0 that an unbiased forecast requires -- so the
'well-calibrated, no systematic over- or under-prediction' reading holds up when it is actually measured.
That matters more than RMSE for risk use: a volatility forecast biased high or low mis-sizes every position
built on it, and the bias would be invisible in the error metric alone.

Neural Networks III — Recurrent Networks & LSTMs¶

The recurrence from scratch → a PyTorch LSTM, vs HAR-RV and GARCH on realized volatility¶

1. The data and the task — realized volatility¶

2. The recurrence — from scratch¶

3. The LSTM — gated memory¶

4. The horse race — LSTM vs the classical volatility models¶

5. Reading the race — the honest verdict¶

6. Proportions vs predictions¶

7. Summary¶