gpin_model.py

# -*- coding: utf-8 -*-

# numpy for matrix algebra
from pandas import Series
from pandas import isnull
import numpy as np
import scipy as sp
from numpy import log, exp, log1p
from math import log as mlog

from scipy.special import gamma, logsumexp, gammaln
import scipy.optimize as op

# import common functions
from common import *

class GPINModel(object):

    def __init__(self,a,r,p,eta,th,d,n=1,t=252):
        """Initializes parameters of EPIN model

        n : the number of stocks to simulate, default 1
        t : the number of periods to simulate, default 252 (one trading year)
        """

        # Assign model parameters
        self.a, self.r, self.p, self.eta, self.th, self.d, self.N, self.T = a, r, p, eta, th, d, n, t

        self.states = self._draw_states()
        self.lamb = _lam(r, p, size=(n, t))

        self.buys = np.random.poisson(self.lamb*((th*(self.states != 1))+((th+eta)*(self.states == 1))))
        self.sells = np.random.poisson(self.lamb*(((1-th)*(self.states != -1))+(((1-th)+eta)*(self.states == -1))))
        self.alpha = compute_alpha(a, r, p, eta, d, th, self.buys, self.sells)

    def _draw_states(self):
        """Draws the states for N stocks and T periods.

        In the Easley and O'Hara sequential trade model at the beginning of each period nature determines whether there is an information event with probability $\alpha$ (a). If there is information, nature determines whether the signal is good news with probability $\delta$ (d) or bad news $1-\delta$ (1-d).

        A quick way to implement this is to draw all of the event states at once as an `NxT` matrix from a binomial distribution with $p=\alpha$, and independently draw all of the news states as an `NxT` matrix from a binomial with $p=\delta$.

        An information event occurs for stock i on day t if `events[i][t]=1`, and zero otherwise. The news is good if `news[i][t]=1` and bad if `news[i][t]=-1`.

        The element-wise product of `events` with `news` gives a complete description of the states for the sequential trade model, where the state variable can take the values (-1,0,1) for bad news, no news, and good news respectively.

        self : EOSequentialTradeModel instance which contains parameter definitions
        """
        events = np.random.binomial(1, self.a, (self.N,self.T))
        news = np.random.binomial(1, self.d, (self.N,self.T))
        news[news == 0] = -1

        states = events*news

        return states

def _lam(r,p,size=None):
    "Compute \lambda_{NI} from shape r and scale p/(1-p) params"
    return np.nan_to_num(np.random.gamma(r,p/(1-p),size))

def _lf(th_b, th_s, r, p, n_buys, n_sells, pdenom=1):
    res =  log(th_b)*n_buys+log(1-th_s)*n_sells - lfact(n_buys) - lfact(n_sells) - gammaln(r) + log(1-p)*r + log(p)*(n_buys+n_sells) + gammaln(r+n_buys+n_sells) - log(pdenom)*r - log(pdenom)*(n_buys+n_sells)
    return res

def _ll(a, r, p, eta, d, th, n_buys, n_sells):
    return np.array([log(1-a)+_lf(th, th, r, p, n_buys, n_sells),
                   log(a*d)+_lf(th+eta, th, r, p, n_buys, n_sells, 1+eta*p),
                   log(a*(1-d))+_lf(th, th-eta, r, p, n_buys, n_sells, 1+eta*p)])

def compute_alpha(a, r, p, eta, d, th, n_buys, n_sells):
    '''Compute the conditional alpha given parameters, buys, and sells.

    '''
    ys = _ll(a, r, p, eta, d, th, n_buys, n_sells)

    ymax = ys.max(axis=0)
    lik = exp(ys-ymax)
    alpha = lik[1:].sum(axis=0)/lik.sum(axis=0)

    return alpha

def nbm_ll(theta, x):
    a,p,eta,r = theta
    q = (p+eta*p)/(1+eta*p)
    def _nbl(a,p,r,x):
        x = x.reset_index(drop=True)
        return log(a)+log(1-p)*r+log(p)*x-lfact(x)-gammaln(r)+gammaln(r+x)
    ll = np.array([_nbl(1-a,p,r,x),_nbl(a,q,r,x)])
    return sum(logsumexp(ll,axis=0))

def _loglik(theta, a, r, p, n_buys, n_sells):
    eta,d,th = theta

    ll = np.array([log(1-a)+_lf(th, th, r, p, n_buys, n_sells),
                   log(a*d)+_lf(th+eta, th, r, p, n_buys, n_sells, 1+eta*p),
                   log(a*(1-d))+_lf(th, th-eta, r, p, n_buys, n_sells, 1+eta*p)])

    return sum(logsumexp(ll,axis=0))


# --- Vectorized likelihood functions for optimizer ---

def loglik(theta, n_buys, n_sells):
    """Vectorized log-likelihood for the full 6-parameter GPIN model.

    Precomputes data-dependent constants and factors out shared terms
    across the 3 states (no-event, good-news, bad-news).
    """
    a, p, eta, r, d, th = theta

    # Data-dependent constants (would be even faster precomputed outside,
    # but scipy.optimize calls this many times with same data)
    total = n_buys + n_sells
    lf_buys = lfact(n_buys)
    lf_sells = lfact(n_sells)

    # Common terms across all 3 states
    gr = gammaln(r)
    log_1mp = mlog(1 - p)
    log_p = mlog(p)
    base = -lf_buys - lf_sells - gr + log_1mp * r + log_p * total

    # State 0: no event — th_b=th, th_s=th, pdenom=1
    log_th = mlog(th)
    log_1mth = mlog(1 - th)
    s0 = mlog(1 - a) + log_th * n_buys + log_1mth * n_sells + base + gammaln(r + total)

    # States 1,2: event — pdenom = 1 + eta*p
    pdenom = 1 + eta * p
    log_pd = mlog(pdenom)
    base_e = base - log_pd * r - log_pd * total + gammaln(r + total)

    # State 1: good news — th_b=th+eta, th_s=th
    s1 = mlog(a * d) + mlog(th + eta) * n_buys + log_1mth * n_sells + base_e

    # State 2: bad news — th_b=th, th_s=th-eta
    s2 = mlog(a * (1 - d)) + log_th * n_buys + mlog(1 - th + eta) * n_sells + base_e

    # logaddexp chain: log(exp(s0) + exp(s1) + exp(s2)) summed over days
    return np.logaddexp(s0, np.logaddexp(s1, s2)).sum()


def loglik_precomp(theta, n_buys, n_sells, lf_buys, lf_sells, total, gr_total):
    """Vectorized log-likelihood with precomputed data constants.

    This is the hot path — called hundreds of times by scipy.optimize.
    All data-dependent terms are precomputed once in fit().
    """
    a, p, eta, r, d, th = theta

    # Parameter-dependent terms
    gr = gammaln(r)
    log_1mp_r = mlog(1 - p) * r
    log_p = mlog(p)

    # Base: terms shared across all 3 states
    # -lfact(B) - lfact(S) + log(1-p)*r + log(p)*(B+S) - gammaln(r) + gammaln(r+B+S)
    base = -lf_buys - lf_sells + log_1mp_r + log_p * total - gr + gammaln(r + total)

    # State 0: no event
    log_th = mlog(th)
    log_1mth = mlog(1 - th)
    s0 = mlog(1 - a) + log_th * n_buys + log_1mth * n_sells + base

    # States 1,2: event with pdenom adjustment
    pdenom = 1 + eta * p
    log_pd = mlog(pdenom)
    adj = -log_pd * (r + total)  # combines -log(pd)*r - log(pd)*total

    # State 1: good news
    s1 = mlog(a * d) + mlog(th + eta) * n_buys + log_1mth * n_sells + base + adj

    # State 2: bad news
    s2 = mlog(a * (1 - d)) + log_th * n_buys + mlog(1 - th + eta) * n_sells + base + adj

    return np.logaddexp(s0, np.logaddexp(s1, s2)).sum()


def nbm_ll_precomp(theta, turn, lf_turn):
    """Vectorized negative binomial mixture log-likelihood with precomputed lfact.

    Two-component NegBin mixture on total turnover.
    """
    a, p, eta, r = theta
    q = (p + eta * p) / (1 + eta * p)

    gr = gammaln(r)
    gr_turn = gammaln(r + turn)
    common = -lf_turn - gr + gr_turn

    # Component 1: no-event intensity
    s0 = mlog(1 - a) + mlog(1 - p) * r + mlog(p) * turn + common

    # Component 2: event intensity
    s1 = mlog(a) + mlog(1 - q) * r + mlog(q) * turn + common

    return np.logaddexp(s0, s1).sum()


def fit(n_buys, n_sells, starts=10, maxiter=100,
        a=None, r=None, p=None, eta=None, th=None, d=None,
        se=None, winsorize_turn=False, **kwargs):
    import pandas as pd
    from statsmodels.regression.linear_model import OLS

    # Convert to numpy once — avoid pandas overhead in optimizer
    n_buys_arr = np.asarray(n_buys, dtype=np.float64)
    n_sells_arr = np.asarray(n_sells, dtype=np.float64)

    turn = n_buys_arr + n_sells_arr
    if winsorize_turn:
        sp.stats.mstats.winsorize(turn, limits=0.05, inplace=True)
    ETA_MAX = (np.abs(n_buys_arr - n_sells_arr) / turn).max()

    # Precompute data-dependent constants (called once, used in every loglik eval)
    lf_buys = lfact(n_buys_arr)
    lf_sells = lfact(n_sells_arr)
    total = n_buys_arr + n_sells_arr
    lf_turn = lfact(turn)

    # Bounds and ranges
    bounds = [(0.00001,0.99999)]*2+[(0.00001,ETA_MAX)]+[(0.00001,np.inf)]+[(0.00001,0.99999)]*2
    ranges = [(0.00001,0.99999)]*2+[(0.00001,ETA_MAX)]+[(0.00001,999)]+[(0.00001,0.99999)]*2

    # Initial parameter estimates
    a0 = a or 0.5
    eta0 = eta or (np.abs(n_buys_arr - n_sells_arr) / turn).mean()
    d0 = d or 0.5
    p0 = p or max(0.00001, 1 - turn.mean() / max(turn.var(), 1e-10))
    r0 = r or max(0.00001, (1 - p0) / p0 * turn.mean())

    results = OLS(n_sells_arr, n_buys_arr).fit()
    th0 = th or (1 / (1 + np.atleast_1d(results.params)[0]))

    # --- Stage 1: NegBin mixture on turnover (4 params) ---
    nll_nbm = lambda theta, *args: -nbm_ll_precomp(theta, *args)
    f_nbm = nll_nbm([a0, p0, eta0, r0], turn, lf_turn)

    for i in range(starts):
        if i:  # random start after first
            x0 = [np.random.uniform(l, np.nan_to_num(h)) for (l, h) in ranges[:4]]
        else:
            x0 = [a0, p0, eta0, r0]
        try:
            res = op.minimize(nll_nbm, x0, method=None,
                              bounds=bounds[:4], args=(turn, lf_turn))
            if res.success and res.fun < f_nbm:
                # Check not on boundary
                on_bound = any(x in b for x, b in zip(res.x, bounds[:4]))
                if not on_bound:
                    f_nbm = res.fun
                    a0, p0, eta0, r0 = res.x
        except Exception:
            continue

    # --- Stage 2: Full 6-param loglik ---
    nll_full = lambda theta, *args: -loglik_precomp(theta, *args)
    f_full = nll_full([a0, p0, eta0, r0, d0, th0],
                      n_buys_arr, n_sells_arr, lf_buys, lf_sells, total, None)
    best_params = [a0, p0, eta0, r0, d0, th0]
    best_hess_inv = None

    for i in range(starts):
        if i:
            x0 = [np.random.uniform(l, np.nan_to_num(h)) for (l, h) in ranges]
        else:
            x0 = [a0, p0, eta0, r0, d0, th0]
        try:
            res = op.minimize(nll_full, x0, method=None,
                              bounds=bounds,
                              args=(n_buys_arr, n_sells_arr, lf_buys, lf_sells, total, None))
            if res.success and res.fun < f_full:
                on_bound = any(x in b for x, b in zip(res.x, bounds))
                if not on_bound:
                    f_full = res.fun
                    best_params = res.x.tolist()
                    best_hess_inv = res.get('hess_inv', None)
        except Exception:
            continue

    a0, p0, eta0, r0, d0, th0 = best_params
    rc = 0

    # --- Output ---
    param_names = 'a,p,eta,r,d,th'.split(',')
    output = dict(zip(param_names + ['f', 'rc'],
                    best_params + [-f_full, rc]))
    if se and best_hess_inv is not None:
        try:
            stderr = 1 / np.sqrt(np.diag(np.linalg.inv(
                best_hess_inv.todense() if hasattr(best_hess_inv, 'todense') else best_hess_inv)))
        except Exception:
            stderr = np.zeros(6)
        output = {'params': dict(zip(param_names, best_params)),
                  'se': dict(zip(param_names, stderr.tolist())),
                  'stats': {'f': -f_full, 'rc': rc}
                 }
    return output

if __name__ == '__main__':

    import pandas as pd
    from regressions import *

    a,r,p,eta,d,th = 0.54, 12, 0.9998, 0.0064, 0.465, 0.469

    N = 1000
    T = 252

    model = GPINModel(a,r,p,eta,th,d,n=N,t=T)

    buys = to_series(model.buys)
    sells = to_series(model.sells)

    aoib = abs(buys-sells)
    poib = abs(sm.OLS(buys, sells, missing='drop').fit().resid)

    turn = buys+sells
    alpha = to_series(model.alpha)

    def run_regs(df):
        # run regression
        m = []
        m.append(partial_r2(df['alpha'],df[['poib','poib2']], df[['poib','poib2','turn','turn2']]))
        out = pd.DataFrame(m, columns=['results'])
        out.index.names = ['model']
        return out

    regtab = pd.DataFrame({'alpha':alpha,'poib':poib,'poib2':poib**2,'turn':turn,'turn2':turn**2})

    res = run_regs(regtab)

    print(est_tab(res.results, est=['params','tvalues'], stats=['rsquared','rsquared_sp']))