The restriction is written as a linear (matrix equation):

$$R_{(m \times k)} B_{(k \times 1)} = q_{(m x 1)}.$$

The test statistic is given by

$$W = (R\hat{\beta}-q)^t(\sigma^2R(X^t X)^{-1}R^t) ^{-1} (R \hat{\beta}-q)$$

where $$W$$ has a $$\chi^2(m)$$ distribution with m degrees of freedom and

standard table of that distribution can be used.

The Null hypothesis $$H_0: R\beta = q$$ is tested against the alternative hypothesis $$H_1: R\beta \ne q$$.

Our python function

`waldtest(R, q, fit)`simply requires the user to input the R and q vectors as single dimensional Python arrays.The fit is the ols fitting object returned by the lm.ols function.

INSERT PYTHON FUNCTION HERE with example!!!

""" """ File waldtest.py Author Dr. Ernesto Adorio U.P. Clark ernesto.adorio@gmail.com description Computes the Wald statistic. version 0.0.1 2010.06.20 first version 0.0.2 2010.06.21 fixed bugs for example. """ from matlib import * from lm import * from dataframe import * def waldtest(R, q, fit): """ Performs a wald test on the linear equality restriction equation RB = q. R is a matrix of dimension m x k. B is a vector of length k. q is a vector of length m. """ m = len(R) B = fit.betas invxtx = fit.VC sigmasqr = fit.MRSS # changed from fit.RSS, 06.21 # compute wald statistic! m, k = matdim(R) if m != len(q): raise ValueError, "waldtest: R must have same rows as length of q!" if k != len(B): raise ValueError,"waldtest: R must have same columns as length of B!" print "debug: sigmasqr=", sigmasqr S = matinverse(matkmul(matABAt(R, invxtx), sigmasqr)) # modified 06.21 RB = matvec(R,B) RBmq = [rb - e for (rb,e) in zip(RB, q)] print "debug:", RBmq W = matvectmatvec(RBmq, S) return W, m def Test(): """ Danao's example page """ data = """ YEAR OUTPUT LABOR CHEM MACH 1947 58 297 15 54 1948 63 285 16 62 1949 62 285 18 68 1950 61 265 19 72 1951 63 251 21 77 1952 66 237 23 81 1953 66 220 24 82 1954 66 214 24 82 1955 69 220 26 83 1956 69 212 27 84 1957 67 196 27 83 1958 73 182 28 83 1959 74 183 32 84 1960 76 177 32 83 1961 76 167 35 80 1962 77 163 38 80 1963 80 155 43 79 1964 79 148 46 80 1965 82 144 49 80 1966 79 132 56 82 1967 83 128 66 85 1968 85 124 69 86 1969 85 118 73 86 1970 84 112 75 85 1971 92 108 81 87 1972 91 110 86 86 1973 93 109 90 90 1974 88 109 92 92 1975 95 106 83 96 1976 97 100 96 98 1977 100 100 100 100 1978 104 100 107 194 1979 111 99 123 104 1980 104 96 123 101 1981 118 96 129 98 1982 116 93 118 94 1983 96 97 105 90 1984 112 98 121 88 1985 113 84 123 83 """ D = dataframe(data) X = matSelectCols(D.data, [2,3,4]) matprint(X) Xlog = [] for x in X: Xlog.append([log(v) for v in x]) matInsertConstCol(Xlog, 0, 1.0) matprint(Xlog) Y = [log(y[1]) for y in D.data] fit = ols(Xlog, Y) # print coefficients. print "Coefficients:", fit.betas R = [[0.,1.,1.,1.]] q = [1] W, m = waldtest(R, q, fit) print "W statistic:", W if __name__ == "__main__": Test()

When the original 0.0.1 program was run, the output was

Coefficients: [2.4262377786012621, 0.10743484613669807, 0.33519197626861796, 0.029383835786802592] W statistic: 0.57599916366

The coefficients match those of Danao, pp. 210, all right. But the W statistic .576 does not match the published value of 4.722226 So we will be in debug mode until we get this right.

After a good sleep, we disccovered that

sigmasqr should be fit.MRSS and the expresion S should be inverted! After these corrections,

the program test outputs the expected 4.722226!

SAFE TO USE THE PYTHON CODE CODE for computing the Wald statistic. but please, ABSOLUTELY NO WARRANTIES.

Will revise the code shortly to return the p-value of the test. Or readers who know Python can help.

We shall use the scipy function to do it.

The file dataframe.py may be found in http://adorio-research.org/wordpress/?p=7834.

We will post the updated matlib.py library soon!

Why does w-stats is one figure hence it is not corresponding the coefficients??

ReplyDeleteTsitso