1"""
2Diagnostics for SUR and 3SLS estimation
3"""
4
5__author__= "Luc Anselin lanselin@gmail.com,    \
6             Pedro V. Amaral pedrovma@gmail.com  \
7             Tony Aburaad taburaad@uchicago.edu"
8
9
10import numpy as np
11import scipy.stats as stats
12import numpy.linalg as la
13from .sur_utils import sur_dict2mat,sur_mat2dict,sur_corr,spdot
14from .regimes import buildR1var,wald_test
15
16
17__all__ = ['sur_setp','sur_lrtest','sur_lmtest','lam_setp','surLMe','surLMlag']
18
19
20def sur_setp(bigB,varb):
21    '''
22    Utility to compute standard error, t and p-value
23
24    Parameters
25    ----------
26    bigB    : dictionary
27              of regression coefficient estimates,
28              one vector by equation
29    varb    : array
30              variance-covariance matrix of coefficients
31
32    Returns
33    -------
34    surinfdict : dictionary
35                 with standard error, t-value, and
36                 p-value array, one for each equation
37
38    '''
39    vvb = varb.diagonal()
40    n_eq = len(bigB.keys())
41    bigK = np.zeros((n_eq,1),dtype=np.int_)
42    for r in range(n_eq):
43        bigK[r] = bigB[r].shape[0]
44    b = sur_dict2mat(bigB)
45    se = np.sqrt(vvb)
46    se.resize(len(se),1)
47    t = np.divide(b,se)
48    tp = stats.norm.sf(abs(t))*2
49    surinf = np.hstack((se,t,tp))
50    surinfdict = sur_mat2dict(surinf,bigK)
51    return surinfdict
52
53def lam_setp(lam,vm):
54    """
55    Standard errors, t-test and p-value for lambda in SUR Error ML
56
57    Parameters
58    ----------
59    lam        : array
60                 n_eq x 1 array with ML estimates for spatial error
61                 autoregressive coefficient
62    vm         : array
63                 n_eq x n_eq subset of variance-covariance matrix for
64                 lambda and Sigma in SUR Error ML
65                 (needs to be subset from full vm)
66
67    Returns
68    -------
69               : tuple
70                 with arrays for standard error, t-value and p-value
71                 (each element in the tuple is an n_eq x 1 array)
72
73    """
74    vvb = vm.diagonal()
75    se = np.sqrt(vvb)
76    se.resize(len(se),1)
77    t = np.divide(lam,se)
78    tp = stats.norm.sf(abs(t))*2
79    return (se,t,tp)
80
81def sur_lrtest(n,n_eq,ldetS0,ldetS1):
82    '''
83    Likelihood Ratio test on off-diagonal elements of Sigma
84
85    Parameters
86    ----------
87    n        : int
88               cross-sectional dimension (number of observations for an equation)
89    n_eq     : int
90               number of equations
91    ldetS0   : float
92               log determinant of Sigma for OLS case
93    ldetS1   : float
94               log determinant of Sigma for SUR case (should be iterated)
95
96    Returns
97    -------
98    (lrtest,M,pvalue) : tuple
99                        with value of test statistic (lrtest),
100                        degrees of freedom (M, as an integer)
101                        p-value
102
103    '''
104    M = n_eq * (n_eq - 1)/2.0
105    lrtest = n * (ldetS0 - ldetS1)
106    pvalue = stats.chi2.sf(lrtest,M)
107    return (lrtest,int(M),pvalue)
108
109
110def sur_lmtest(n,n_eq,sig):
111    '''
112    Lagrange Multiplier test on off-diagonal elements of Sigma
113
114    Parameters
115    ----------
116    n        : int
117               cross-sectional dimension (number of observations for an equation)
118    n_eq     : int
119               number of equations
120    sig      : array
121               inter-equation covariance matrix for null model (OLS)
122
123    Returns
124    -------
125    (lmtest,M,pvalue) : tuple
126                        with value of test statistic (lmtest),
127                        degrees of freedom (M, as an integer)
128                        p-value
129    '''
130    R = sur_corr(sig)
131    tr = np.trace(np.dot(R.T,R))
132    M = n_eq * (n_eq - 1)/2.0
133    lmtest = (n/2.0) * (tr - n_eq)
134    pvalue = stats.chi2.sf(lmtest,M)
135    return (lmtest,int(M),pvalue)
136
137
138def surLMe(n_eq,WS,bigE,sig):
139    """
140    Lagrange Multiplier test on error spatial autocorrelation in SUR
141
142    Parameters
143    ----------
144    n_eq       : int
145                 number of equations
146    WS         : array
147                 spatial weights matrix in sparse form
148    bigE       : array
149                 n x n_eq matrix of residuals by equation
150    sig        : array
151                 cross-equation error covariance matrix
152
153    Returns
154    -------
155    (LMe,n_eq,pvalue) : tuple
156                        with value of statistic (LMe), degrees
157                        of freedom (n_eq) and p-value
158
159    """
160    # spatially lagged residuals
161    WbigE = WS * bigE
162    # score
163    EWE = np.dot(bigE.T,WbigE)
164    sigi = la.inv(sig)
165    SEWE = sigi * EWE
166    #score = SEWE.sum(axis=1)
167    #score.resize(n_eq,1)
168    # note score is column sum of Sig_i * E'WE, a 1 by n_eq row vector
169    # previously stored as column
170    score = SEWE.sum(axis=0)
171    score.resize(1,n_eq)
172
173
174    # trace terms
175    WW = WS * WS
176    trWW = np.sum(WW.diagonal())
177    WTW = WS.T * WS
178    trWtW = np.sum(WTW.diagonal())
179    # denominator
180    SiS = sigi * sig
181    Tii = trWW * np.identity(n_eq)
182    tSiS = trWtW * SiS
183    denom = Tii + tSiS
184    idenom = la.inv(denom)
185    # test statistic
186    #LMe = np.dot(np.dot(score.T,idenom),score)[0][0]
187    # score is now row vector
188    LMe = np.dot(np.dot(score,idenom),score.T)[0][0]
189    pvalue = stats.chi2.sf(LMe,n_eq)
190    return (LMe,n_eq,pvalue)
191
192def surLMlag(n_eq,WS,bigy,bigX,bigE,bigYP,sig,varb):
193    """
194    Lagrange Multiplier test on lag spatial autocorrelation in SUR
195
196    Parameters
197    ----------
198    n_eq       : int
199                 number of equations
200    WS         : spatial weights matrix in sparse form
201    bigy       : dictionary
202                 with y values
203    bigX       : dictionary
204                 with X values
205    bigE       : array
206                 n x n_eq matrix of residuals by equation
207    bigYP      : array
208                 n x n_eq matrix of predicted values by equation
209    sig        : array
210                 cross-equation error covariance matrix
211    varb       : array
212                 variance-covariance matrix for b coefficients (inverse of Ibb)
213
214    Returns
215    -------
216    (LMlag,n_eq,pvalue) : tuple
217                          with value of statistic (LMlag), degrees
218                          of freedom (n_eq) and p-value
219
220    """
221    # Score
222    Y = np.hstack((bigy[r]) for r in range(n_eq))
223    WY = WS * Y
224    EWY = np.dot(bigE.T,WY)
225    sigi = la.inv(sig)
226    SEWE = sigi * EWY
227    score = SEWE.sum(axis=0)  # column sums
228    score.resize(1,n_eq)  # score as a row vector
229
230    # I(rho,rho) as partitioned inverse, eq 72
231    # trace terms
232    WW = WS * WS
233    trWW = np.sum(WW.diagonal()) #T1
234    WTW = WS.T * WS
235    trWtW = np.sum(WTW.diagonal()) #T2
236
237    # I(rho,rho)
238    SiS = sigi * sig
239    Tii = trWW * np.identity(n_eq) #T1It
240    tSiS = trWtW * SiS
241    firstHalf = Tii + tSiS
242    WbigYP = WS * bigYP
243    inner = np.dot(WbigYP.T, WbigYP)
244    secondHalf= sigi * inner
245    Ipp= firstHalf + secondHalf #eq. 75
246
247    # I(b,b) inverse is varb
248
249    # I(b,rho)
250    bp = sigi[0,] * spdot(bigX[0].T,WbigYP)  #initialize
251    for r in range(1,n_eq):
252        bpwork = sigi[r,] * spdot(bigX[r].T,WbigYP)
253        bp = np.vstack((bp,bpwork))
254    # partitioned part
255    i_inner = Ipp - np.dot(np.dot(bp.T, varb), bp)
256    # partitioned inverse of information matrix
257    Ippi = la.inv(i_inner)
258
259    # test statistic
260    LMlag = np.dot(np.dot(score,Ippi),score.T)[0][0]
261    # p-value
262    pvalue = stats.chi2.sf(LMlag,n_eq)
263    return (LMlag,n_eq,pvalue)
264
265def sur_chow(n_eq,bigK,bSUR,varb):
266    """
267    test on constancy of regression coefficients across equations in
268    a SUR specification
269
270    Note: requires a previous check on constancy of number of coefficients
271    across equations; no other checks are carried out, so it is possible
272    that the results are meaningless if the variables are not listed in
273    the same order in each equation.
274
275    Parameters
276    ----------
277    n_eq       : int
278                 number of equations
279    bigK       : array
280                 with the number of variables by equation (includes constant)
281    bSUR       : dictionary
282                 with the SUR regression coefficients by equation
283    varb       : array
284                 the variance-covariance matrix for the SUR regression
285                 coefficients
286
287    Returns
288    -------
289    test       : array
290                 a list with for each coefficient (in order) a tuple with the
291                 value of the test statistic, the degrees of freedom, and the
292                 p-value
293
294    """
295    kr = bigK[0][0]
296    test = []
297    bb = sur_dict2mat(bSUR)
298    kf = 0
299    nr = n_eq
300    df = n_eq - 1
301    for i in range(kr):
302        Ri = buildR1var(i,kr,kf,0,nr)
303        tt,p = wald_test(bb,Ri,np.zeros((df,1)),varb)
304        test.append((tt,df,p))
305    return test
306
307def sur_joinrho(n_eq,bigK,bSUR,varb):
308    """
309    Test on joint significance of spatial autoregressive coefficient in SUR
310
311    Parameters
312    ----------
313    n_eq       : int
314                 number of equations
315    bigK       : array
316                 n_eq x 1 array with number of variables by equation
317                 (includes constant term, exogenous and endogeneous and
318                 spatial lag)
319    bSUR       : dictionary
320                 with regression coefficients by equation, with
321                 the spatial autoregressive term as last
322    varb       : array
323                 variance-covariance matrix for regression coefficients
324
325    Returns
326    -------
327               : tuple
328                 with test statistic, degrees of freedom, p-value
329
330    """
331    bb = sur_dict2mat(bSUR)
332    R = np.zeros((n_eq,varb.shape[0]))
333    q = np.zeros((n_eq,1))
334    kc = -1
335    for i in range(n_eq):
336        kc = kc + bigK[i]
337        R[i,kc] = 1
338    w,p = wald_test(bb,R,q,varb)
339    return (w,n_eq,p)
340