Solving Ax = b s.t. x contains all positive values

Specific Domains Numerics
,
1

I am solving an asset pricing equations via discretizing an underlying Markov process, X \in R^d, resulting in \tilde{X} \sim \Pi, where \Pi is a transition matrix whose rows sum to unity, and all entries in [0,1].

Let \delta \in R^d be a constant vector s.t. \eta_t := \exp \left(\delta X_t \right) represents growth rate of discounted cash flows. Let \tilde{\eta_i} := \exp\left(\delta \tilde{X}_i \right) be its discretized counterpart.

The asset pricing equation in the continuous space is:

V(x) = \eta_t + \mathbb{E}[\eta_t \eta_{t+1} | X_t = x] + \mathbb{E}[\eta_t \eta_{t+1} \eta_{t+2}| X_t = x] + \cdots

Cleary V(x) >0 \ \forall x \in R^d, because \eta >0, and the solution is a stochastic integral.

The above asset-pricing equation has a recursive structure:

V(x) = \eta_t + \eta_t \mathbb{E}[V(X_{t+1}) | X_t = x].

For the discretized process \tilde{X}, the recursive structure results in a system of linear equations,

\vec{V} = \vec{\eta} + diag(\vec{\eta}) \Pi \vec{V},

whose solution is

\vec{V} = \left( \mathbb{I} - diag\left(\vec{\eta}\right) \Pi\right)^{-1} \vec{\eta}.

The issue I’m having is using the backslash operator to solve the systems of equations sometimes results in a solution where \vec{V} has large negative entries, even though domain knowledge tells me all entries should be positive.

Is there a way in Julia to solve a system of linear equations s.t. all elements of the solution are greater than 0?

Also, any tips on solving the indefinite stochastic integral are welcome. In my application X is a Markov-switching vector autoregression where D = 4.

2

You can use JuMP: Introduction · JuMP

The exact formulation depends on whether you expect there to be an exact solution for Ax = b, or whether by A \ b you’re really finding some minimum norm solution.

One approach is:

function solve(A::Matrix, b::Vector)
    m, n = size(A)
    using JuMP, HiGHS
    model = Model(HiGHS.Optimizer)
    @variable(model, x[1:n] >= 0)
    @constraint(model, A * x .== b)
    optimize!(model)
    return value.(x)
end

Another approach is:

function solve(A::Matrix, b::Vector)
    m, n = size(A)
    using JuMP, HiGHS
    model = Model(HiGHS.Optimizer)
    @variable(model, x[1:n] >= 0)
    @variable(model, residuals[1:m])
    @constraint(model, residuals .== A * x .- b)
    @objective(model, Min, sum(residuals.^2))
    optimize!(model)
    return value.(x)
end
3

Non-negative least squares would probably help too. NNLS.jl or NonNegLeastSquares.jl would be good places to start

1 Like
4

Thank you, both.

I found the pivot algorithm in NonNegLeastSquares.jl to be most performant on my type fo problem (Large and sparse).

1 Like
5

I think you have an analysis problem here first, not an implementation problem.

Is that system singular (in exact arithmetic)? Otherwise it will have a unique solution; is it always positive in theory? Can you prove it? If the answer is no, NNLS will put a band-aid on it and return something positive, but it will not be an exact solution to the system and to your problem.

Your description is a little fuzzy so I cannot make out the details: what are X_t, \tilde{X}_i and how are they related to \Pi? Should there be a \tilde{\eta} in the discretized version (last two equations) instead of \eta? If you fix it, we might be able to tell you more about the exact solution to that linear system and figure out if there is a problem in the theory or in the implementation.