Solving linear system in Julia: Mismatch with solution found using R

Question

I need to find a vector z such that Az = b;

That is, z = inverse(A) x b.

A = a 5 x 5 matrix.

B = a 5 x 1 vector.

I had earlier solved (correctly) for z in R. But now I am using Julia and can’t seem to get the correct answer.

R code:

# a 3 x 3 variance-covariance matrix of asset returns

s <- matrix(c(0.0100, 0.0018, 0.0011, 0.0018, 0.0109, 0.0026, 0.0011, 0.0026, 0.0199), nrow=3, ncol=3)

# a vector of 3 mean returns

m =  c(0.0427, 0.0015, 0.0285)

## Create the matrix A using S and m ## as given in the theory.

top = cbind(2*s, m, rep(1, 3))  ## creates top 3 rows of  matrix A
mid = c(m, 0, 0)            ## creates 4th row of matrix A
bot = c(rep(1, 3), 0, 0)        ## creates 5th row of matrix A
A = rbind(top, mid, bot)

A

##   Output:                             
    0.0200  0.0036      0.0022      0.0427      1
    0.0036  0.0218      0.0052      0.0015      1
    0.0022  0.0052      0.0398      0.0285      1
    0.0427  0.0015      0.0285      0.0000      0
    1.0000  1.0000      1.0000      0.0000      0

# Create the vector b using the first asset’s mean

b = c(rep(0, 3), m[1], 1)

b
## Output:  0.0000      0.0000      0.0000      0.0427      1.0000

z = solve(A)%*%b

z

## Outputs: 
   0.827454555
  -0.090746123
   0.263291567
 -0.350290278
  -0.001844252

All of these are in fact the correct outputs, as verified from another source.

But when I replicate the same process in Julia I cannot get the same solution for z !!

Julia code:

# vector of ones (to be used in A below)
n_assets = 3
o = ones(Int8, n_assets)

Data for means and std.devs.
rA = .0427 ;    sA = .1000;     lA = "A" # For asset A etc.
rB = .0015 ;     sB = .1044;     lB = "B"
rC = .0285 ;     sC = .1411;     lC = "C"

varA = sA^2;    varB = sB^2;    varC = sC^2;
covAB = .0018;      covBC = .0026;      covAC = .0011
covBA = covAB;      covCB = covBC;      covCA = covAC

# vector of means
μ = [rA rB rC] ## a row vector       # same as m in the R code
# create var-covar matrix
Σ = [varA covAB covAC; covBA varB covBC; covCA covCB varC]  # same as s in the R code

### find A matrix

A_top = hcat(2*Σ, μ', o) ## o is a col vector of ones
A_mid = hcat(μ,0,0 )
A_bot = hcat(o',0,0)
A = vcat(A_top, A_mid, A_bot)

μp0 = rA            ## .0427 from the data
b0 = vcat(z, μp0, 0)        # z is a col vector of zeros

## solve for z in Az = b0
z = inv(A)*b0
z
## outputs:
5-element Vector{Float64}:
  0.9263700069697688
 -1.0942794778158993
  0.16790947084613073
 -0.8399181116023101
  0.020907108510299026

Why am I not getting the correct z vector that I am getting in R? I tried rounding off the var-cover matrix to 4 decimal digits but to no effect. I checked the var-covar matrix, its inverse, A, its inverse and the b vectors and they all are matching in both cases. It is just the z vector that is not correct in Julia output. Am I missing something elementary here?

You didn't print `A` in Julia. The most likely kind of error in a problem like this is loading the matrix incorrectly; e.g. maybe it ended up transposed compared to R. Another minor comment: rather than `solve(A) %*% b`, it's better to use the more efficient `solve(A, b)` in R. I imagine there's something equivalent in Julia. — user2554330, Sep 08 '21 at 19:56

Peace Wang · Accepted Answer · 2021-09-09T01:54:31.163

Your b0 in julia is different with its value in R. In Julia, the last number of b0 is 0, which should be 1 (b0 = vcat(zeros(3), μp0, 1)).

b0 = vcat(zeros(3), μp0, 0)

 0.0
 0.0
 0.0
 0.0427
 0.0

A straightforard check is

A = [0.0200  0.0036      0.0022      0.0427      1;
    0.0036  0.0218      0.0052      0.0015      1;
    0.0022  0.0052      0.0398      0.0285      1;
    0.0427  0.0015      0.0285      0.0000      0;
    1.0000  1.0000      1.0000      0.0000      0]

b = [0.0000;      0.0000 ;     0.0000;      0.0427;      1.0000]

inv(A)*b
# or 
# A/b is faster

5-element Vector{Float64}:
  0.8274545552796227
 -0.09074612277886507
  0.2632915674992422
 -0.35029027763746184
 -0.001844251656967253

Solving linear system in Julia: Mismatch with solution found using R

1 Answers1