Open in app

Sign in

Write

Sign in

Analysis of a Naive Statistical Arbitrage for Trading

Adam
Call For Atlas
Published in
8 min readDec 3, 2023
Basket of stocks done by MidJourney 12.2023

Statistical arbitrage (StatArb) is any technique in quantitative finance that uses statistical and mathematical models to exploit a short-term market inefficiency. A known StatArb is Pairs trading, which we discussed in a previous article and utilizes cointegration to trade a net market-neutral portfolio of two assets.

In this version of the StatArb algorithm, we will use a basket of weighted stocks (leading stocks and a target) to capture temporary lead-lag relationships in the market.

Prepare your Environment

Have a jupyter environment ready, and pip install these libraries:

  • numpy
  • pandas
  • yfinance

You’ll need access to analysis_utils library for common functions.

Set up the Universe

We will load a basket of the highest capitalized stocks from the Nasdaq 100.

Since Stat Arb is a tactical strategy, we will focus on correlation and covariance only, and filter the universe of stocks by this criteria.

START_DATE = "2020-01-01"
END_DATE = "2023-11-22"
tickers = ["AAPL", "MSFT", "AMZN", "TSLA", "GOOG", "GOOGL", "META", "NVDA", "PEP", "AVGO", "ADBE", "COST", "PYPL", "AMD", "QCOM", "INTC", "TXN", "CHTR", "TMUS", "ISRG", "SBUX", "AMGN", "INTU", "ADP", "CSX", "ADI", "MU", "ZM", "MAR", "GILD", "MELI", "WDAY", "CTSH", "PANW", "REGN", "LRCX", "BKNG", "EXC", "ADSK", "KLAC", "TEAM"]

tickers_df = load_ticker_prices_ts_df(tickers, START_DATE, END_DATE)
tickers_rets_df = tickers_df.dropna(axis=1).pct_change().dropna()  # first % is NaN
# 1+ to allow the cumulative product of returns over time, and -1 to remove it at the end.
tickers_rets_df = (1 + tickers_rets_df).cumprod() - 1
plt.figure(figsize=(16, 16))
for ticker in tickers_rets_df.columns:
    plt.plot(tickers_rets_df.index, tickers_rets_df[ticker] * 100.0, label=ticker)
plt.xlabel("Date (Year-Month)")
plt.ylabel("Cummulative Returns(%")
plt.legend()
plt.show()
TARGET = "AAPL"
MA_WINDOW = 24
ARB_WINDOW = MA_WINDOW * 10
plt.figure(figsize=(16, 8))

corr_ticks = []
LEADS = []
for ticker in tickers:
    if ticker == TARGET:
        continue
    correlation = tickers_df[ticker].corr(tickers_df[TARGET])
    if abs(correlation) < 0.75:
        continue
    LEADS.append(ticker)
    corr_ts = tickers_df[ticker].rolling(ARB_WINDOW).corr(tickers_df[TARGET])
    plt.plot(
        tickers_df.index,
        corr_ts.rolling(MA_WINDOW).mean(),
        label=f"{ticker} (Corr: {correlation:.2f})",
        alpha=correlation,
        linewidth=correlation,
    )
    corr_ticks.append(ticker)
plt.axhline(y=0, color="k", linestyle="--", linewidth=1)
plt.xlabel("Date (Year-Month)")
plt.ylabel("Correlation with AAPL")
plt.legend()
plt.show()

The above plot looks at the correlations overtime, we see that these oscillate between uncorrelated and correlated.

In actual implementations, the basket of stocks gets redone periodically on specific factors — correlation could be one of these. Luckily for us, most of the chosen tickers like MSFT and GOOGL remain highly correlated during this simulation.

Signalling Action

Let’s go through the pseudocode for this naive StatArb algo we are building:

  1. First we compute the exponential moving average (EMA) using a set window (we test 1 month). This is done for all stocks both leads and the target.
  2. We calculate the delta of the current price from the EMA for all stocks.
  3. Using a rolling window (we test 1/2 a year), we get the pair-wise correlation between the EMA differences of the leads and target stock. This captures the pairs short-term trends.
  4. We do the same for the covariance ratios between each pair. This quantifies how strong the leader’s price changes can project the target’s.
  5. We use the covariance to make projections. Each leader’s own EMA difference is multiplied by its covariance ratio with the target to project what the target’s difference “should” be, if it moved in lockstep with the leader. This is similar to the cointegration algo in pairs-trading.
  6. The error between the projected and actual target EMA difference is calculated. This residual indicates how much the target is under or over performing relative to the leader-based prediction. The sum of these will be our signal.
  7. Finally we use the correlation ratios to weight each lead’s input to the signal.

Just as a reminder, covariance is:

  • Xbar and X is the mean of X and an individual observation.
  • Ybar and Y is the mean of Y and an individual observation.
  • N is the number of observations

We use the covariance to give us the magnitude of the movement of the leading stock, not the correlation coefficients (although this normalized value has some magnitude in it, this information has been lost). By itself covariance has no inherent interpretation.

here is the code:

arb_df = tickers_rets_df.copy()
arb_df["price"] = tickers_df[TARGET]
arb_df[f"ema_{TARGET}"] = arb_df[TARGET].ewm(MA_WINDOW).mean()
arb_df[f"ema_d_{TARGET}"] = arb_df[TARGET] - arb_df[f"ema_{TARGET}"]
arb_df[f"ema_d_{TARGET}"].fillna(0, inplace=True)

for ticker in LEADS:
    arb_df[f"ema_{ticker}"] = arb_df[ticker].ewm(MA_WINDOW).mean()
    arb_df[f"ema_d_{ticker}"] = arb_df[ticker] - arb_df[f"ema_{ticker}"]
    arb_df[f"ema_d_{ticker}"].fillna(0, inplace=True)
    arb_df[f"{ticker}_corr"] = (
        arb_df[f"ema_d_{ticker}"].rolling(ARB_WINDOW).corr(arb_df[f"ema_d_{TARGET}"])
    )
    arb_df[f"{ticker}_covr"] = (
        arb_df[[f"ema_d_{ticker}", f"ema_d_{TARGET}"]]
        .rolling(ARB_WINDOW)
        .cov(numeric_only=True)
        .groupby(level=0, axis=0, dropna=True)  # Cov returns pairwise!
        .apply(lambda x: x.iloc[0, 1] / x.iloc[0, 0])
    )
    arb_df[f"{ticker}_emas_d_prj"] = (
        arb_df[f"ema_d_{ticker}"] * arb_df[f"{ticker}_covr"]
    )
    arb_df[f"{ticker}_emas_act"] = (
        arb_df[f"{ticker}_emas_d_prj"] - arb_df[f"ema_d_{TARGET}"]
    )
arb_df.filter(regex=f"(_emas_d_prj|_corr|_covr)$").dropna().iloc[ARB_WINDOW:]

We now weigh the signals:

ts = 0
delta_projected = 0
for ticker in LEADS:
    corr_abs = abs(arb_df[f"{ticker}_corr"].fillna(0))
    weights += corr_abs
    arb_df[f"{ticker}_emas_act_w"] = arb_df[f"{ticker}_emas_act"].fillna(0) * corr_abs
    delta_projected += arb_df[f"{ticker}_emas_act_w"]

weights = weights.replace(0, 1)
arb_df[f"{TARGET}_signal"] = delta_projected / weights
arb_df[f"{TARGET}_signal"].iloc[ARB_WINDOW:]

To visualize the quality of our signal, we’ll calculate the errors and the signals confidence intervals at 95%.

Where:

  • Xhat is the sample mean,
  • Z is the Z-score corresponding to the desired confidence level (1.96 for a 95% confidence interval),
  • sigma is the populations’ standard deviation or in our case, the residuals.
  • n is the sample size.

The code below will calculate these criteria, and will be used later when acting on the signal:

errors = (
    arb_df[f"ema_{TARGET}"] + arb_df[f"ema_d_{TARGET}"] + arb_df[f"{TARGET}_signal"]
) - (arb_df[f"ema_{TARGET}"].shift(-1))
arb_df["rmse"] = np.sqrt((errors**2).rolling(ARB_WINDOW).mean()).fillna(0)

me = errors.rolling(ARB_WINDOW).mean().fillna(0)
e_std = errors.rolling(ARB_WINDOW).std().fillna(0)
ci = (me - 1.96 * e_std, me + 1.96 * e_std)
arb_df["ci_lower"] = arb_df[f"{TARGET}_signal"].fillna(0) + ci[0]
arb_df["ci_upper"] = arb_df[f"{TARGET}_signal"].fillna(0) + ci[1]
arb_df["ci_spread"] = arb_df["ci_upper"] - arb_df["ci_lower"]
arb_df.fillna(0, inplace=True)
arb_df[["ci_lower", "ci_upper"]].iloc[ARB_WINDOW:].tail(10)

Visualizing the Signal and Errors

With the signal and error metrics ready, let’s visualize it all:

fig, axes = plt.subplots(2, 1, gridspec_kw={"height_ratios": (3, 1)}, figsize=(22, 8))

axes[0].set_title(f"Visualizing {TARGET}'s Signal and Errors")
axes[0].plot(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df[TARGET].iloc[ARB_WINDOW:],
    label=f"{TARGET} $",
    alpha=1,
    color="b",
)
axes[0].plot(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df[f"ema_{TARGET}"].iloc[ARB_WINDOW:],
    label=f"{TARGET} EMA $",
    alpha=1,
)
axes[0].plot(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df[f"{TARGET}_signal"].iloc[ARB_WINDOW:].fillna(0)
    + arb_df[TARGET].iloc[ARB_WINDOW:],
    label=f"{TARGET} + Signal $",
    alpha=0.75,
    linestyle="--",
    color="r",
)
axes[0].legend()
axes[1].plot(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df[f"{TARGET}_signal"].iloc[ARB_WINDOW:],
    label=f"Wieghted {TARGET} Signal $",
    alpha=0.75,
    color="g",
)
axes[1].plot(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df["rmse"].iloc[ARB_WINDOW:],
    label="RMSE",
    alpha=0.75,
    linestyle="--",
    color="r",
)
axes[1].fill_between(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df["ci_lower"].iloc[ARB_WINDOW:],
    arb_df["ci_upper"].iloc[ARB_WINDOW:],
    color="gray",
    alpha=0.3,
)
axes[1].axhline(0, color="black", linestyle="--", linewidth=1)
axes[1].legend()
plt.tight_layout()
plt.show()

As the StatArb’s window rolls across more data, it’s predictions become more ‘accurate’. The tighter the confidence bands, the more likely we will simulate a trade.

I do have to say that in the post-covid market regime (2020–2023), there were some long running trends, helping our predictions’ accuracy.

Simulated Trading

Now let’s act on our signals. This is a simulation, therefore we are ignoring various market nuances, including stop losses and transaction expenses.

the code below will long/short a portfolio depending on the signal’s position to selected thresholds:

LONG_THRESHOLD = 0.0025
SHORT_THRESHOLD = -0.0025
CONF_SPREAD_THRESHOLD = 0.15  
MAX_SHARES = 1

arb_df["orders"] = 0
signals = arb_df[f"{TARGET}_signal"]
prev_signals = signals.shift(-1)
add_long_cond = (signals > LONG_THRESHOLD) & (prev_signals <= LONG_THRESHOLD) & (signals < arb_df["ci_upper"]) & (signals > arb_df["ci_lower"]) & (arb_df["ci_spread"] < CONF_SPREAD_THRESHOLD)
add_short_cond = (signals < SHORT_THRESHOLD) & (prev_signals >= SHORT_THRESHOLD) & (signals < arb_df["ci_upper"]) & (signals > arb_df["ci_lower"]) & (arb_df["ci_spread"] < CONF_SPREAD_THRESHOLD)

arb_df.loc[add_long_cond, "orders"] += MAX_SHARES
arb_df.loc[add_short_cond, "orders"] -= MAX_SHARES
arb_df["orders"].fillna(0, inplace=True)
arb_df.loc[arb_df["orders"] != 0, "orders"].tail(10)

This code section uses Pandas to derive a PnL. A real life trade-execution system would be much different than this:

signal_changes_df = arb_df.loc[(add_long_cond | add_short_cond), ["price", "orders", f"{TARGET}_signal"]]

signal_changes_df["holdings"] = signal_changes_df["orders"].cumsum()
signal_changes_df["stat_chng"] = np.sign(signal_changes_df["orders"].shift(1)) != np.sign(signal_changes_df["orders"])
prev_holdings = signal_changes_df["holdings"].shift(1)
signal_changes_df["price_open"] = signal_changes_df["price"].shift(1)
signal_changes_df["cost_open_avg"] = (signal_changes_df.loc[signal_changes_df["stat_chng"] == False, "price_open"].shift(1) + signal_changes_df.loc[signal_changes_df["stat_chng"]== False, "price_open"]) / 2
signal_changes_df["cost_open_avg"].fillna(signal_changes_df["price_open"], inplace=True)
signal_changes_df["price_close"] = signal_changes_df["price"]
signal_changes_df["price_close"].iloc[0] = np.nan  # First signal shouldn't have a closing price
signal_changes_df["pnl"] = (signal_changes_df.loc[signal_changes_df["stat_chng"], "price_close"]  - signal_changes_df.loc[signal_changes_df["stat_chng"], "cost_open_avg"]) * np.sign(signal_changes_df["holdings"].shift(1))
signal_changes_df["pnl_rets"] = signal_changes_df["pnl"] / signal_changes_df["cost_open_avg"].abs()
signal_changes_df.fillna(0, inplace=True)
arb_df = pd.concat([arb_df, signal_changes_df[["pnl", "pnl_rets", "holdings"]]], axis=1).drop_duplicates(keep='last')
signal_changes_df.tail(10)

Finally we visualize the trade’s executed with the signal and the acceptable error measures:

plt.figure(figsize=(26, 18))

fig, (ax1, ax2, ax3) = plt.subplots(
    3, 1, gridspec_kw={"height_ratios": (3, 1, 1)}, figsize=(24, 12)
)
ax1.plot(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df["price"].iloc[ARB_WINDOW:],
    color="g",
    lw=1.25,
    label=f"{TARGET}",
)
ax1.plot(
    arb_df.loc[add_long_cond].index,
    arb_df.loc[add_long_cond, "price"],
    "^",
    markersize=12,
    color="blue",
    label="Buy",
)
ax1.plot(
    arb_df.loc[add_short_cond].index,
    arb_df.loc[add_short_cond, "price"],
    "v",
    markersize=12,
    color="red",
    label="Sell",
)
ax1.set_ylabel("Price in $")
ax1.legend(loc="upper left", fontsize=10)
ax1.set_title(f"Stat-Arb {TARGET}", fontsize=18)
ax2.plot(
    arb_df["pnl"].iloc[ARB_WINDOW:].index, arb_df["pnl_rets"].iloc[ARB_WINDOW:].fillna(0).cumsum(), color="b", label="Returns"
)
ax2.set_ylabel("Cumulative Profit (%)")
ax2.legend(loc="upper left", fontsize=10)
ax2.set_title(f"Cummulative Returns", fontsize=18)
ax3.plot(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df[f"{TARGET}_signal"].iloc[ARB_WINDOW:],
    label=f"Signal",
    alpha=0.75,
    color="g",
)
ax3.plot(
    arb_df.loc[add_long_cond].index,
    arb_df.loc[add_long_cond, f"{TARGET}_signal"],
    "^",
    markersize=12,
    color="blue",
    label="Buy",
)
ax3.plot(
    arb_df.loc[add_short_cond].index,
    arb_df.loc[add_short_cond, f"{TARGET}_signal"],
    "v",
    markersize=12,
    color="red",
    label="Sell",
)
ax3.axhline(0, color="black", linestyle="--", linewidth=1)
ax3.axhline(SHORT_THRESHOLD, color="r", linestyle="--", label=f"Short ({SHORT_THRESHOLD})",linewidth=1)
ax3.axhline(LONG_THRESHOLD, color="b", linestyle="--", label=f"Long ({LONG_THRESHOLD})",linewidth=1)
ax3.fill_between(arb_df.index, arb_df[f"{TARGET}_signal"], SHORT_THRESHOLD, where=(arb_df[f"{TARGET}_signal"] < SHORT_THRESHOLD), interpolate=True, color='red', alpha=0.3)
ax3.fill_between(arb_df.index, arb_df[f"{TARGET}_signal"], LONG_THRESHOLD, where=(arb_df[f"{TARGET}_signal"] > LONG_THRESHOLD), interpolate=True, color='blue', alpha=0.3)
ax3.fill_between(
    arb_df.iloc[ARB_WINDOW:].index,
    arb_df["ci_lower"].iloc[ARB_WINDOW:],
    arb_df["ci_upper"].iloc[ARB_WINDOW:],
    color="gray",
    alpha=0.3,
    label=f"Confidence"
)
ax3.legend(loc="lower left", fontsize=10)
ax3.set_title(f"Leads' Signal", fontsize=18)
plt.tight_layout()
plt.show()

As this is a simulation, we have simplified the market as we did the StatArb Algo itself. A real world variant would use sophisticated statistical processes in a high-frequency execution environment, trading intraday or close to this, as alpha from market inefficiencies is taken in an instant.

We also skipped the back-testing, forward-testing, and parameter tuning of parameters like the rolling window lengths, and thresholds.

Conclusion

In this article, we simulated a naive statistical arbitrage trading strategy, built around confidence intervals and variance-weighted signals.

As most finfluencing stuff out there, this is a toy strategy and should not be used beyond this analysis. Stay safe.

lock stepping with the market by MidJourney 12.2023

References

Github

Article here is also available on Github

Kaggle notebook available here

Media

All media used (in the form of code or images) are either solely owned by me, acquired through licensing, or part of the Public Domain and granted use through Creative Commons License.

CC Licensing and Use

This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

Algorithmic Trading
Statistical Arbitrage
Quantitative Finance
Python
Machine Learning

--

--

Adam
Call For Atlas

Written by Adam

480 Followers
Editor for

People, tech, and product focused senior technology leader with a penchant for AI in quant finance: https://www.linkedin.com/in/adam-darmanin/ #CallForAtlas

Help

Status

About

Careers

Press

Blog

Privacy

Terms

Text to speech

Teams