Re: st: RE: Matrix question- square roots of matrix elements

 From [email protected] (Richard Gates) To [email protected] Subject Re: st: RE: Matrix question- square roots of matrix elements Date Thu, 25 Jan 2007 06:52:37 -0600

```To add to Kit Baum's statalist post, Mata has LAPACK under the hood so another
approach to the least squares problem is to use QR decomposition:

clear
qui {
input y x1 x2 x3
123.1   1       1.92    12.4
124.3   1       2.15    9.9
89.3    1       1.67    2.4
141.3   1       1.68    13.8
112.8   1       1.75    3.5
108.1   1       1.55    1.8
143.9   1       1.54    17.8
124.2   1       2.1     9.8
110.1   1       2.44    8.3
111.7   1       2.47    9.8
123.8   1       1.86    12.6
123.5   1       1.93    11.5
110.2   1       2.47    7.4
100.9   1       2.11    6.1
123.3   1       2.1     9.5
115.7   1       1.73    8.8
116.6   1       1.86    4.9
153.5   1       2.19    18.8
149.2   1       1.9     18.9
89      1       1.67    2.3
132.6   1       2.43    14.1
97.5    1       2.13    2.9
106.1   1       2.33    5.9
115.3   1       1.75    7.6
98.5    1       2.05    5.3
135.1   1       2.35    16.8
124.2   1       2.12    8.8
98.4    1       2.13    3.2
114.8   1       1.89    5.4
142.5   1       1.5     17.3
122.6   1       1.93    11.2
127.7   1       2.27    11.2
113     1       1.66    7.9
144.2   1       1.73    17
109.2   1       1.59    3.3
106.8   1       2.29    7.1
145     1       1.86    15.3
124     1       1.91    12.7
106.7   1       2.34    6.1
153.2   1       2.13    19.6
120.1   1       2.05    6.3
119.3   1       1.89    9
150.6   1       2.12    18.7
92.2    1       1.87    2.2
130.5   1       2.09    16
112.5   1       1.76    4.5
111.8   1       1.77    4.3
120.1   1       1.94    9.3
107.4   1       2.37    8.3
128.6   1       2.1     15.4
124.6   1       2.29    9.2
127.2   1       2.36    10.2
end
}
mata

st_view(y=.,.,"y")
st_view(X=.,.,("x1", "x2" ,"x3"))

H = tau = R = J(0,0,0)
p = (1,0,0)

hqrdp(X,H,tau,R,p)
qy = hqrdmultq(H,tau,y,1)

r = rank(R)
R1 = R[|1,1\r,r|]

p1 = invorder(p[1::r])

b = solveupper(R,qy[1::r])[p1]
Rinv = solveupper(R,I(r))
sse = sum(qy[(r+1)::length(qy)]:^2)
ssr = sum(qy[2::r]:^2)
s2 = sse/(length(qy)-r)
se = sqrt(s2:*diagonal(Rinv*Rinv'))[p1]
t = b:/se

r2 =  ssr/(sse+ssr)

(b,se,t)
sqrt(s2)
(ssr,sse)
r2
end

regress y x2 x3

Note that the estimated condition number of X'X is considerablly larger than
that of X by itself.  The condition number is a measure of the linear
system sensitivity.  The smaller the number the better.

mata
: cond(X)
94.62867081

: cond(cross(X,X))
8954.585338

: m = colsum(X[,2::3]):/rows(X)
: cond(crossdev(X[,2::3],m,X[,2::3],m))
369.2486364

-Rich
rgates@stata
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```