Skip to content

Commit

Permalink
Qr (#253)
Browse files Browse the repository at this point in the history 
* remove cholesky

* speed up

* update

* rmv find_cols

* patched version

* Update Project.toml

* Update fit.jl
  • Loading branch information
@matthieugomez
matthieugomez authored Nov 30, 2023
1 parent 9e5d020 commit 79a372b
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 35 deletions.
24 changes: 13 additions & 11 deletions src/fit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,16 +321,17 @@ function StatsAPI.fit(::Type{FixedEffectModel},
perm = 1:(n_exo + n_endo)

# first pass: remove colinear variables in Xendo
notcollinear_fe_endo = collinear_fe[find_cols_endo(n_exo, n_endo)] .== false
basis_endo = basis(eachcol(Xendo)...; has_intercept = false) .* notcollinear_fe_endo
notcollinear_fe_endo = collinear_fe[(n_exo+2):(n_exo+n_endo+1)] .== false
basis_endo = basis(Xendo; has_intercept = false) .* notcollinear_fe_endo
Xendo = getcols(Xendo, basis_endo)

# second pass: remove colinear variable in Xexo, Z, and Xendo
notcollinear_fe_exo = collinear_fe[find_cols_exo(n_exo)] .== false
notcollinear_fe_z = collinear_fe[find_cols_z(n_exo, n_endo, n_z)] .== false
notcollinear_fe_exo = collinear_fe[2:(n_exo+1)] .== false
notcollinear_fe_z = collinear_fe[(n_exo+n_endo+2):(n_exo+n_endo+n_z+1)] .== false
notcollinear_fe_endo_small = notcollinear_fe_endo[basis_endo]

basis_all = basis(Xexo, Z, eachcol(Xendo)...; has_intercept = has_intercept)

basis_all = basis(Xexo, Z, Xendo; has_intercept = has_intercept)
basis_Xexo = basis_all[1:size(Xexo, 2)] .* notcollinear_fe_exo
basis_Z = basis_all[(size(Xexo, 2) +1):(size(Xexo, 2) + size(Z, 2))] .* notcollinear_fe_z
basis_endo_small = basis_all[(size(Xexo, 2) + size(Z, 2) + 1):end] .* notcollinear_fe_endo_small
Expand Down Expand Up @@ -380,7 +381,7 @@ function StatsAPI.fit(::Type{FixedEffectModel},
# get linearly independent columns
n_exo = size(Xexo, 2)
perm = 1:n_exo
notcollinear_fe_exo = collinear_fe[find_cols_exo(n_exo)] .== false
notcollinear_fe_exo = collinear_fe[2:(n_exo+1)] .== false
basis_Xexo = basis(Xexo; has_intercept = has_intercept) .* notcollinear_fe_exo
Xexo = getcols(Xexo, basis_Xexo)
Xhat = Xexo
Expand All @@ -393,6 +394,7 @@ function StatsAPI.fit(::Type{FixedEffectModel},
## Do the regression
##
##############################################################################

crossx = Xhat' * Xhat
Xy = Symmetric(hvcat(2, crossx, Xhat'y, zeros(size(Xhat, 2))', [0.0]))
invsym!(Xy; diagonal = 1:size(Xhat, 2))
Expand Down Expand Up @@ -466,15 +468,15 @@ function StatsAPI.fit(::Type{FixedEffectModel},
# Compute Fstat of First Stage
if has_iv && first_stage
Pip = Pi[(size(Pi, 1) - size(Z_res, 2) + 1):end, :]
#try
try
r_kp = Vcov.ranktest!(Xendo_res, Z_res, Pip,
vcov_method, size(X, 2), dof_fes)
p_kp = chisqccdf(size(Z_res, 2) - size(Xendo_res, 2) + 1, r_kp)
F_kp = r_kp / size(Z_res, 2)
# catch
# @info "ranktest failed ; first-stage statistics not estimated"
# p_kp, F_kp = NaN, NaN
# end
catch
@info "ranktest failed ; first-stage statistics not estimated"
p_kp, F_kp = NaN, NaN
end
end

# Compute rss, tss
Expand Down
35 changes: 11 additions & 24 deletions src/utils/basecol.jl
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,15 @@
# create the matrix X'X
function crossprod(xs...)
idx = vcat(0, cumsum(size(x, 2) for x in xs))
XX = zeros(idx[end], idx[end])
for i in 1:length(xs)
for j in i:length(xs)
XX[(idx[i]+1):idx[i+1], (idx[j]+1):idx[j+1]] .= xs[i]' * xs[j]
end
end
return Symmetric(XX, :U)
crossprod(x1) = Symmetric(x1'x1)

function crossprod(x1, x2)
Symmetric(hvcat(2, x1'x1, x1'x2,
zeros(size(x2, 2), size(x1, 2)), x2'x2))
end

function crossprod(x1, x2, x3)
Symmetric(hvcat(3, x1'x1, x1'x2, x1'x3,
zeros(size(x2, 2), size(x1, 2)), x2'x2, x2'x3,
zeros(size(x3, 2), size(x1, 2)), zeros(size(x3, 2), size(x2, 2)), x3'x3))
end

# generalized 2inverse
Expand Down Expand Up @@ -47,7 +49,6 @@ function getcols(X::AbstractMatrix, basecolX::AbstractVector)
sum(basecolX) == size(X, 2) ? X : X[:, basecolX]
end


function ls_solve(X, y::AbstractVector)
Xy = crossprod(X, y)
invsym!(Xy, diagonal = 1:size(X, 2))
Expand All @@ -59,17 +60,3 @@ function ls_solve(X, Y::AbstractMatrix)
invsym!(XY, diagonal = 1:size(X, 2))
return XY[1:size(X, 2),(end-size(Y, 2)+1):end]
end

##############################################################################
# Auxiliary functions to find columns of exogeneous, endogenous and IV variables
##############################################################################

function find_cols_exo(n_exo)
2:n_exo+1
end
function find_cols_endo(n_exo, n_endo)
n_exo+2:n_exo+n_endo+1
end
function find_cols_z(n_exo, n_endo, n_z)
n_exo+n_endo+2:n_exo+n_endo+n_z+1
end

2 comments on commit 79a372b

@matthieugomez
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/96258

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.10.0 -m "<description of version>" 79a372bcf41e60be7f1680bfc7d79ef6dd942eb0
git push origin v1.10.0

Please sign in to comment.