<commit_before>"""
Maximum likelihood estimation of spatial process models


TODO
----

 - check ORD versus FULL likelihood results with sparse implementation for ORD
 - asymptotic standard errors for ML_Lag
 - ML_error likelihood
 - ML_error partial
 - ML_error estimation

"""

__author__ = "<NAME> <EMAIL>, <NAME> <EMAIL>"

import numpy as np
import pysal as ps
import scipy.sparse as SPARSE


def defl_lag(r, w, e1, e2):
    """
    Derivative of the likelihood for the lag model

    Parameters
    ----------

    r: estimate of rho

    w: spatial weights object

    e1: ols residuals of y on X 

    e2: ols residuals of Wy on X

    Returns
    -------

    tuple: 
        tuple[0]: dfl (scalar) value of derivative at parameter values
        tuple[1]: tr (scalar) trace of  (I-rW)^{-1} W

    """
    n = w.n
    e1re2 = e1 - r * e2
    num = np.dot(e1re2.T, e2)
    den = np.dot(e1re2.T, e1re2)
    a = -r * w.full()[0]
    np.fill_diagonal(a, np.ones((w.n,1)))
    a = np.linalg.inv(a)
    tr = np.trace(np.dot(a, w.full()[0]))
    dfl = n * (num/den) - tr
    return dfl[0][0], tr

def log_like_lag_full(w,b,X,y):
    n = w.n
    r = b[0]    # ml estimate of rho
    b = b[1:]   # ml for betas
    ldet = _logJacobian(w, r)
    yl = ps.lag_spatial(w,y)
    ys = y - r * yl
    XX = np.dot(X.T, X)
    iXX = np.linalg.inv(XX)
    b = np.dot(iXX, np.dot(X.T,ys))
    yhat = r * yl + np.dot(X,b)
    e = y - yhat
    e2 = (e**2).sum()
    sig2 = e2 / n
    ln2pi = np.log(2*np.pi)
    return ldet - n/2. * ln2pi - n/2. * np.log(sig2) - e2/(2 * sig2)


def c_log_like_lag_full(r, e1, e2, w):
    """
    Log concentrated likelihood for the lag model using full evaluation
    """
    n = w.n
    e1re2 = e1 - r*e2
    ldet = _logJacobian(w, r)

    return -(n/2.)  * np.log( (e1re2**2).sum() / n  ) + ldet

def log_like_lag_ord(w, b, X, y, evals):
    n = w.n
    r = b[0]    # ml estimate of rho
    b = b[1:]   # ml for betas
    yl = ps.lag_spatial(w,y)
    ys = y - r * yl
    XX = np.dot(X.T, X)
    iXX = np.linalg.inv(XX)
    b = np.dot(iXX, np.dot(X.T,ys))
    yhat = r * yl + np.dot(X,b)
    e = y - yhat
    e2 = (e**2).sum()
    sig2 = e2 / n
    ln2pi = np.log(2*np.pi)
    revals = r * evals
    ldet = np.log(1-revals).sum()
    return ldet - n/2. * ln2pi - n/2. * np.log(sig2) - e2/(2 * sig2)



def c_log_like_lag_ord(r, e1, e2, evals):
    """
    Log concentrated likelihood for the lag model using Ord (1975)
    approximation

    Parameters
    ----------

    r: estimate of rho

    e1: ols residuals of y on X 

    e2: ols residuals of Wy on X

    evals: vector of eigenvalues of W


    Returns<commit_msg>- centralizing log likelihood for the lag model<commit_after>"""
Maximum likelihood estimation of spatial process models


TODO
----

 - check ORD versus FULL likelihood results with sparse implementation for ORD
 - asymptotic standard errors for ML_Lag
 - ML_error likelihood
 - ML_error partial
 - ML_error estimation

"""

__author__ = "<NAME> <EMAIL>, <NAME> <EMAIL>"

import numpy as np
import pysal as ps
import scipy.sparse as SPARSE


def defl_lag(r, w, e1, e2):
    """
    Derivative of the likelihood for the lag model

    Parameters
    ----------

    r: estimate of rho

    w: spatial weights object

    e1: ols residuals of y on X 

    e2: ols residuals of Wy on X

    Returns
    -------

    tuple: 
        tuple[0]: dfl (scalar) value of derivative at parameter values
        tuple[1]: tr (scalar) trace of  (I-rW)^{-1} W

    """
    n = w.n
    e1re2 = e1 - r * e2
    num = np.dot(e1re2.T, e2)
    den = np.dot(e1re2.T, e1re2)
    a = -r * w.full()[0]
    np.fill_diagonal(a, np.ones((w.n,1)))
    a = np.linalg.inv(a)
    tr = np.trace(np.dot(a, w.full()[0]))
    dfl = n * (num/den) - tr
    return dfl[0][0], tr

def log_like_lag_full(w,b,X,y):
    r = b[0]
    ldet = _logJacobian(w, r)
    return log_like_lag(ldet, w, b, X, y)


def c_log_like_lag_full(r, e1, e2, w):
    """
    Log concentrated likelihood for the lag model using full evaluation
    """
    n = w.n
    e1re2 = e1 - r*e2
    ldet = _logJacobian(w, r)

    return -(n/2.)  * np.log( (e1re2**2).sum() / n  ) + ldet

def log_like_lag_ord(w, b, X, y, evals):
    r = b[0]
    revals = r * evals
    ldet = np.log(1-revals).sum()
    return log_like_lag(ldet, w, b, X, y)

def log_like_lag(ldet, w, b, X, y):
    n = w.n
    r = b[0]    # ml estimate of rho
    b = b[1:]   # ml for betas
    yl = ps.lag_spatial(w,y)
    ys = y - r * yl
    XX = np.dot(X.T, X)
    iXX = np.linalg.inv(XX)
    b = np.dot(iXX, np.dot(X.T,ys))
    yhat = r * yl + np.dot(X,b)
    e = y - yhat
    e2 = (e**2).sum()
    sig2 = e2 / n
    ln2pi = np.log(2*np.pi)
    return ldet - n/2. * ln2pi - n/2. * np.log(sig2) - e2/(2 * sig2)



def c_log_like_lag_ord(r, e1, e2, evals):
    """
    Log concentrated likelihood for the lag model using Ord (1975)
    approximation

    Parameters
    ----------

    r: estimate of rho

    e1: ols residuals of y on X 

    e2: ols residuals of Wy on X

    evals: vector of eigenvalues of W


    Returns
{'fraction_non_alphanumeric': 0.14700078926598265, 'fraction_numerical': 0.02584846093133386, 'mean_word_length': 2.0908536585365853, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3780226', 'n_tokens_mistral': 2173, 'n_tokens_neox': 2075, 'n_words': 611}