How to compute the combined probability of loss for 2 time series (consisting of historical stock prices)? May I please ask the community's support with the following problem?
I have 2 time series, with approximately 1000 observations each (same number of observations for both). They represent the daily closing prices for 2 stocks: asset A and asset B. Assuming we are currently at time t, I'm interested in computing (forecasting) the probability that, (at anytime) over the next 5 days (i.e. t+1, t+2, t+3, t+4, t+5), either stock A's or B's price will fall by at least $r_{loss}$%:
$$
Pr[P^{A}_{t+i}<P^{A}_{t}*(1-r_{loss})\ \textbf{or}\ P^{B}_{t+i}<P^{B}_{t}*(1-r_{loss})\ for\ some\ i\in\{1,2,3,4,5\}]=\ ?
$$

For added clarity:

*

*t = current time; I only have information (observations) up to
and including time t

*$P^{A}, P^{B}$ = price of stock A, price of stock  B

*$r_{loss}$ = loss threshold, $r_{loss}\in[0,1]$. For example, $r_{loss}=0.03$ represents a $3\%$ loss.


Here is an approach based on the ideas that:

*

*the logs of stock prices can be modeled with Brownian motion

*in Brownian motions $$Pr[\max_{i\le n} X_i>t] = 2 Pr[X_n>t]$$ by the reflection principle (https://en.m.wikipedia.org/wiki/Reflection_principle_(Wiener_process))

*in a standard normal distribution, the probability $p$ of getting a draw below $-3/\sqrt{5}$ is substantial ($\sim9\%$), but $p^2$ is negligible

These suggest the following procedure:

*

*Calculate the series of daily log-returns for each stock, $R_t =\ln(P_t/P_{t-1})$.

*Calculate the standard deviations $\sigma_A$ and $\sigma_B$ of those log-returns, and $\sigma_m=\max(\sigma_A,\sigma_B)$.

*Count the (backward-looking) fraction of days $t$ so far where either $A$ or $B$ went below $P_t e^{-3\sigma_m}$ over the next five days; call this fraction $L_b$.

*Find the $\rho$ for which $$L_b = 2Pr[\min(X,Y)<-3\sigma_m]$$ when $X$ and $Y$ have a bivariate normal distribution with means $0$, variances $5\sigma_A^2$ and $5\sigma_B^2$ and correlation $\rho$. (This probability is given by an integral, and the appropriate $\rho$ can be found by trial and error or more sophisticated root-finding.)

*[Update] Expanding $Pr[\min(X,Y)<t]$ as a quadratic Taylor series in $\rho$ gives the nice form
$$1-\Phi(u)\Phi(v)-\frac{2\rho+\rho^2uv}{4\pi}\exp\left(-\frac{\sigma_A^2+\sigma_B^2}{2\sigma_A\sigma_B}uv\right)$$
where $u=-t/(\sqrt{5}\sigma_A)$, $v=-t/(\sqrt{5}\sigma_B)$, and that may be good enough for this method.

*Let $$L_f = 2Pr[\min(X,Y)<\ln(1-r_{loss})]$$ where $X$ and $Y$ have that same distribution with the same $\rho$.

Then $L_f$ is a prediction for the (forward-looking) fraction of days where either $A$ or $B$ will lose at least $r_{loss}$ over the next five days.
We can check that this gives good answers in several cases:

*

*when $1-r_{loss}=e ^{-3\sigma_m}$, since then $L_f=L_b$

*when $r_{loss}=1$, since then $L_f= 0$

*when $r_{loss}=-\infty$ (i.e. allowing any possible gains), since then $L_f\simeq 1$

*when $\sigma_A$ is much larger than $\sigma_B$, when $A$ and $B$ are independent, or when $A$ and $B$ are perfectly correlated

*if the $5$ is replaced by another number of days or an infinite number of days

There are surely other approaches which can better capture other features of the stock prices, but this may be the simplest solution that works as well as it does.
EDIT: code by mihnea_11235
def prob_of_future_loss(historical_prices,  nr_std = 3, r_loss = .03, period = 5):
    """
    Function that computes L_f
    :param historical_prices:
    :param nr_std:
    :param r_loss:
    :param period:
    :return:
    """

    def get_returns(data):
        """
        Helper function to get log-returns

        :param data: historical prices
        :return:
        """
        returns_all = np.log(data) - np.log(data.shift(1))
        # drop 1st row
        returns_all = returns_all.dropna(how='any')
        return returns_all


    def get_Lb(prices, std_max, nr_std, period):
        """
        Function that computes the (backward-looking) fraction of days t so far where either series_1 or series_2 went
        below P_{t}exp(−3*std_m) over the next five days; call this fraction Lb

        :param prices:
        :param std_max:
        :return:
        """
        asset_1 = prices.columns[0]
        asset_2 = prices.columns[1]

        # variable to count if 1 or 0
        count = []
        for idx in range(0, prices.shape[0]-period):

            # check if P_{t+i}<P_{t}*exp(-nr_std*std_max) for  some i = {1,2,3,...,period}
            check_1 = np.where(prices[asset_1].iloc[idx+1 : idx + period + 1] < prices[asset_1].iloc[idx] * np.exp(-1 * nr_std * std_max))[0]
            check_2 = np.where(prices[asset_2].iloc[idx+1 : idx + period + 1] < prices[asset_2].iloc[idx] * np.exp(-1 * nr_std * std_max))[0]

            if len(check_1)>0 or len(check_2)>0:
                count.append(1)
            else:
                count.append(0)

        return np.sum(count)/len(count)


    def get_Rho(L_b, t, std_1, std_2, period):
        """
        Function that computes the correlation coefficient (rho), given L_b (the (backward-looking) fraction of days
        t so far where either series_1 or series_2 went below P_{t}exp(−3*std_m) over the next five days;

        :param L_b:
        :param t:
        :param std_1:
        :param std_2:
        :param std_max:
        :return:
        """

        # compute interm parts of equations
        u = (-1 * t) / (np.sqrt(period) * std_1)
        v = (-1 * t) / (np.sqrt(period) * std_2)

        # computed cdfs
        cdf_u = norm.cdf(u)
        cdf_v = norm.cdf(v)

        # compute exp part
        E = np.exp(-1 * ((std_1**2 + std_2**2) / (2*std_1*std_2)) * u * v)

        # express as 2n degree equation and find roots
        coef_1 = (E * u * v)/(4*np.pi)
        coef_2 = E / (2 * np.pi)
        coef_3 = cdf_u * cdf_v - 1 + L_b / 2

        # get roots
        coefs = np.roots([coef_1,coef_2,coef_3])

        print('Roots: {}'.format(coefs))

        # use largest root for correlation?
        rho = np.max(coefs)

        return rho


    def get_L_f(rho, t, std_1, std_2, period):
        """
        Function that computes the (forward-looking) fraction of days where either Asset 1 or Asset 2 will lose at
        least r_loss over the next period (e.g. 5) days

        :param rho:
        :param t: np.log(1-r_loss)
        :param std_1:
        :param std_2:
        :return:
        """
        # compute interm parts of equations
        u = -1 * t / (np.sqrt(period) * std_1)
        v = -1 * t / (np.sqrt(period) * std_2)

        # computed cdfs
        cdf_u = norm.cdf(u)
        cdf_v = norm.cdf(v)

        # compute exp part
        E = np.exp(-1 * ((std_1**2 + std_2**2) / (2*std_1*std_2)) * u * v)

        # rho_part
        rho_part = (2 * rho + rho**2 * u * v)/(4*np.pi)

        # compute L_f
        L_f = 2*(1 - cdf_u * cdf_v - rho_part * E)


        return L_f


    # Step 1: get log prices of each series
    returns = get_returns(data=historical_prices)
    asset_1_name = returns.columns[0]
    asset_2_name = returns.columns[1]

    # Step2: compute standard deviations of each series
    std_1 = returns.std()[asset_1_name]#np.std(returns[asset_1_name].values)
    std_2 = returns.std()[asset_2_name]#np.std(returns[asset_2_name].values)
    # get max standard deviation
    std_max = np.max([std_1,std_2])

    # Step 3: find the (backward-looking) fraction of days t so far where either  or  went below P_{t}exp(−3*std_m)
    # over the next five days; call this fraction Lb
    L_b = get_Lb(prices = historical_prices, std_max = std_max, period = period, nr_std =nr_std)
    print('Lb = {}'.format(L_b))

    # Step 4: find rho given L_b
    rho = get_Rho(L_b = L_b, t = -1 * nr_std * std_max, std_1 = std_1, std_2 = std_2, period = period)

    # Step 5: find L_f given rho,
    L_f = get_L_f(rho = rho, t = np.log(1-r_loss), std_1 = std_1, std_2 = std_2, period = period)
    print('Lf = {}'.format(L_f))

    return L_f

EDIT: Code addendum by Matt F.
The above code by mihnea_11235 can be run in Python by installing yfinance, and then using code such as
import yfinance as yf
import numpy as np
from scipy.stats import norm
history = yf.download("SPY GLD", start="2018-12-26", end="2021-12-15")
prob_of_future_loss(history['Adj Close'])

