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/
