Dear Julia Users, I am currently developing/converting several Bayesian spatial econometric models from MATLAB into Julia. I have successfully coded the spatial autoregressive model (SAR) with diffuse priors in Julia but have a question regarding the use of sparse matrix algorithms. The models use a spatial weight matrix which is usually sparse in nature and this matrix appears in many of the expressions, especially in the random walk Metropolis algorithm used to obtain the marginal distribution for the spatial autocorrelation parameter rho. Here is a code snippet and the complete code is attached:
# Draw for rho
A = (In - rho*W);
denominator = log(det(A))-.5*sigma^-2*dot(A*y-x*beta,A*y-x*beta);
accept = 0;
rho_proposed = rho + zta*randn(1);
while accept == 0
if ((rho_proposed[1] > -1) & (rho_proposed[1] < 1));
accept = 1;
else
rho_proposed = rho + zta*randn(1);
end
end
B = (In - rho_proposed[1]*W);
numerator = log(det(B))-.5*sigma^-2*dot(B*y-x*beta,B*y-x*beta);
u = rand(1);
if ((numerator[1] - denominator[1]) > exp(1))
pp = 1;
else
ratio = exp(numerator[1] - denominator[1]);
pp = min(1,ratio);
end
if (u[1] < pp[1])
rho = rho_proposed[1];
acc = acc + 1;
end
ar = acc/gibbs;
if ar < 0.4
zta = zta/1.1;
end
if ar > 0.6
zta = zta*1.1;
end
A = (In - rho*W);
Basically, I'm wondering if there is any way to make the "A" and "B" matrices
sparse and possibly make it run faster, especially in the log(det(A)) terms.
Currently, 110 draws (with 10 burn) takes approximately 8 seconds on my 64 bit,
core i7 laptop. The computational speed decreases with the sample size n
because the weight matrices are treated as full.
Any help would be greatly appreciated and if anyone is interested in running
the code, the Distributions and Distance packages must be included and
initiated first.
Regards,
Don
sar_gibbs.jl
Description: Binary data
knn_weight_matrix.jl
Description: Binary data
