Economics Example

[azev77]

Hi and thank you for this package.

Solving optimal control problems is at the heart of economics.
(including LQR/LQG etc)
I'm working to see what packages in the Julia ecosystem can be helpful.
See my note on discourse.

I want to replicate that (simple) firm investment problem with your package.

Or consider a slightly more general economics HW problem from Ken Kasa:
Continuous Time

Discrete Time

PS: here are my closed form solutions

[azev77]

I tried solving the firm investment problem with your package:

The only difference is my P is your docs Q, and vice versa...

Now in Julia

julia> using MatrixEquations, LinearAlgebra;

julia> δ=0.1; r=0.0; z=0.27;

julia> A=[1.0 0.0; 0.0 (1.0 - δ)];

julia> B=[0.0 0.0; 0.0 1.0];

julia> Q=[0.0 0.0; z 0.0];

julia> R=[0.0 0.0; -1.0 -0.5];

julia> P, CLSEIG, F = ared(A,B,R,Q);
ERROR: DimensionMismatch("R must be a symmetric/hermitian matrix of dimension 2")
Stacktrace:
 [1] gared(A::Matrix{Float64}, E::UniformScaling{Bool}, B::Matrix{Float64}, R::Matrix{Float64}, Q::Matrix{Float64}, S::Matrix{Float64}; as::Bool, rtol::Float64)
   @ MatrixEquations C:\Users\azevelev\.julia\packages\MatrixEquations\mIx3W\src\riccati.jl:790
 [2] #ared#28
   @ C:\Users\azevelev\.julia\packages\MatrixEquations\mIx3W\src\riccati.jl:693 [inlined]
 [3] ared (repeats 2 times)
   @ C:\Users\azevelev\.julia\packages\MatrixEquations\mIx3W\src\riccati.jl:693 [inlined]
 [4] top-level scope
   @ REPL[8]:1

Is it possible to solve this simple Riccati equation with your package?

PS: I just verified by hand that iteration converges

julia> P0=zeros(size(A))
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

julia> P0 = Q +(A')*P0*A -(A')*P0*B*pinv(R + (B')*P0*B)*(B')*P0*A
2×2 Matrix{Float64}:
 0.0   0.0
 0.27  0.0

julia> P0 = Q +(A')*P0*A -(A')*P0*B*pinv(R + (B')*P0*B)*(B')*P0*A
2×2 Matrix{Float64}:
 0.0    0.0
 0.513  0.0

julia> P0 = Q +(A')*P0*A -(A')*P0*B*pinv(R + (B')*P0*B)*(B')*P0*A
2×2 Matrix{Float64}:
 0.0     0.0
 0.7317  0.0

julia> P0 = Q +(A')*P0*A -(A')*P0*B*pinv(R + (B')*P0*B)*(B')*P0*A
2×2 Matrix{Float64}:
 0.0      0.0
 0.92853  0.0

# a bunch more iterations...

julia> P0 = Q +(A')*P0*A -(A')*P0*B*pinv(R + (B')*P0*B)*(B')*P0*A
2×2 Matrix{Float64}:
 0.0      0.0
 2.69996  0.0

julia> P0 = Q +(A')*P0*A -(A')*P0*B*pinv(R + (B')*P0*B)*(B')*P0*A
2×2 Matrix{Float64}:
 0.0      0.0
 2.69997  0.0

[baggepinnen]

You'd probably need to reformulate the problem into a minimization problem, and thus change the sign of R. You might also have reversed R and Q in the call to ared, since the reference you follow appears to use a different convention for the variable name than in control. Lastly, make your R matrix symmetric by

R = (R+R') ./ 2

[andreasvarga]

I think there are some notational confusions.

The only difference is my P is your docs Q, and vice versa...

You probably meant R and Q, since P is the solution of the Riccati equation.

For the solvability of the discrete-time Riccati equation, Q and R must be (at least) symmetric and for stability, the pair (A,B) must be stabilizable (i.e., must exist a state-feedback matrix F such that A-BF has all eigenvalues in the unit circle centered in the origin). Such an F does not exist for your problem, because the eigenvalue of A at 1 cannot be modified by feedback (check A-B*rand(2,2)).
Note that the stabilizing solution P is always symmetric (which is not the case for P0).

Regarding the LQR problem: The weigthed problems can be easily reformulated as standard problems by replacing A with A+ρ/2I in continuous-time or A with sqrt(β)*A and B with sqrt(β)*B
in the discrete-time.

The convergence you reported can just happen (it is not a proof of anything). By the way: How you intend to use the converged value P0?

My final remark: If the LQR problem is properly formulated, I cannot imagine any reason why it can not be solved using Riccati equation based approach (at least not for this size of the problem).

[azev77]

1. change the sign of R (max to min problem) {don't think this changes anything, not sure why?}
2. reverse order of R, Q in ared() {I think I reversed the order correctly}
3. make symmetric R = (R+R') ./ 2. DONE for R, Q
4. reweighted A, B with sqrt(β)*A , sqrt(β)*B. DONE (this was important)

Below I compute the closed form P_sol and compare it to the MatrixEquations.jl solution.
The numerical P has the same value function slope as P_sol, but the wrong value function intercept...
How can this be improved?

My problem w/ revised notation:

PS: here is a note applying LQ to economics.
Here is their Julia code.

using MatrixEquations, LinearAlgebra; 
δ=0.1; r=0.05; z= r + δ +.02; β=1.0/(1.0 + r); # parameters

I_SS = (z-r-δ)/(r+δ);
K_SS = I_SS/δ;
D_SS = z*K_SS - I_SS - 0.5*(I_SS)^(2.0);
V_SS = ((1.0+r)/r)*D_SS;
Q_SS = (1+r)*(1+I_SS);

P_sol = [(V_SS - Q_SS*K_SS) 0.0; Q_SS   0.0]; # closed form P. Value = x'*P*x
F_sol = [-1.0 0.0; -I_SS   0.0];              # closed form u. Policy = -F_sol*x

A=[1.0 0.0; 0.0  (1.0 - δ)];
B=[0.0 0.0; 0.0   1.0];
Q=[0.0 0.0; z       0.0];
R=[0.0 0.0; -1.0 -0.5];

R = -1.0*R;    # not sure how it helps
Q = -1.0*Q;

A = sqrt(β)*A; B = sqrt(β)*B;             # scale by discount factor
Q = (Q+Q') ./ 2.0; R = (R+R') ./ 2.0;  # make symmetric  
P_sol = (P_sol + P_sol') ./ 2.0; 
# F_sol = (F_sol + F_sol') / 2.0          # F should not be symmetric, not quadratic form etc...

                      
P, CLSEIG, F = ared(A,B,R,Q);            # solve for numerical value function = x'*P*x
F_n  = inv(R + (B')*P*B) * (B')*P*A;      # numerical policy = -F_n*x. SAME as F from ared().

julia> P
2×2 Matrix{Float64}:
 1.08579e-15  0.595
 0.595        1.14451e-16

julia> P_sol
2×2 Matrix{Float64}:
 0.186667  0.595
 0.595     0.0

# the off-diagonal term is correct. main diagonal is wrong...


julia> F_n
2×2 Matrix{Float64}:
 -1.13333      -1.96202e-16
  1.88738e-16   3.26742e-32

julia> F_sol
2×2 Matrix{Float64}:
 -1.0       0.0
 -0.133333  0.0

# Policy function wrong...
# Top left term in F_n equal to sum of two left terms in F_sol???


#P = -1.0*P;

[azev77]

Home & away from a computer now.

My best guess: the problem needs to be formulated differently.

Above I had:
Augmented state vector w/ two states
K: true state
1: a dummy state =1

Augmented control vector w/ two controls
I: true control
1: dummy control =1

I believe the problem can be rewritten without a dummy control variable by adding the matrix ’S’ to the objective.

Is there a result in control theory that says I shouldn't create a dummy control vector as I did above?

[andreasvarga]

I had a look to the note applying LQ to economics. Although I had no time to understand the details, but you should be careful to ensure that the conditions listed there are fulfilled (i.e., R and Q must be symmetric, otherwise the solution can not result symmetric). The trick (R+R')/2 can be applied only if R is "almost" symmetric (e.g., R is the result of operations with floating point matrices). In your case, this is certainly not working. You must choose an Q and R which are symmetric (and possibly also possitive definite, albeit not a necessary condition).

[baggepinnen]

The trick (R+R')/2 can be applied only if R is "almost" symmetric (e.g., R is the result of operations with floating point matrices). In your case, this is certainly not working.

Why wouldn't that work? The trick does not change the value of any quadratic form x'Rx:

julia> x = randn(2);

julia> R = [1 0; 1 2]
2×2 Matrix{Int64}:
 1  0
 1  2

julia> x'R*x
5.164061450407917

julia> R = (R + R') ./ 2;

julia> x'R*x
5.164061450407918

[azev77]

@andreasvarga I'm not sure I understand.
The matrixes Q, R, P only appear in quadratic forms.
For example:
Make symmetric @baggepinnen
We get the same quadratic form:

using MatrixEquations, LinearAlgebra; 
δ=0.1; r=0.05; z= r + δ +.02; β=1.0/(1.0 + r); # parameters

I_SS = (z-r-δ)/(r+δ);
K_SS = I_SS/δ;
D_SS = z*K_SS - I_SS - 0.5*(I_SS)^(2.0);
V_SS = ((1.0+r)/r)*D_SS;
Q_SS = (1+r)*(1+I_SS);

P_sol = [(V_SS - Q_SS*K_SS) 0.0; Q_SS   0.0]; # closed form P. Value = x'*P*x
F_sol = [-I_SS   0.0];                        # closed form u. Policy = -F_sol*x

A=[1.0 0.0; 0.0  (1.0 - δ)];
B=[0.0; 1.0];
Q=[0.0 0.0; z 0.0];
R=[-0.5];
S=[-1.0; 0.0];

Q = -1.0*Q;    # max to min
R = -1.0*R;
S = -1.0*S;

A = sqrt(β)*A; # scale by discount factor
B = sqrt(β)*B;             

Q = (Q + Q') ./ 2.0;   # make symmetric  
R = (R + R') ./ 2.0;  
P_sol = (P_sol + P_sol') ./ 2.0; 
                      
P, CLSEIG, F = ared(A,B,R,Q,S);             # solve for numerical value function = x'*P*x
Fn  = pinv(R .+ (B')*P*B) * ((B')*P*A +S'); # numerical policy = -F_n*x. SAME as F from ared().

# Compare Value function P
# the off-diagonal term is correct (except should be positive). main diagonal is wrong...
julia> P
2×2 Matrix{Float64}:
 -7.88667  -0.595
 -0.595    -0.0

julia> P_sol
2×2 Matrix{Float64}:
 0.186667  0.595
 0.595     0.0

# Compare Policy function F. u = -F*x
# Left term in F equal to 1 + left term in F_sol???
julia> F
1×2 Matrix{Float64}:
 0.866667  0.0

julia> F_sol
1×2 Matrix{Float64}:
 -0.133333  0.0

Repeat w/o multiplying Q,R,S by -1.0

# Compare Value function P
# the off-diagonal term is correct. main diagonal is wrong...
julia> P
2×2 Matrix{Float64}:
 7.88667  0.595
 0.595    0.0

julia> P_sol
2×2 Matrix{Float64}:
 0.186667  0.595
 0.595     0.0

# Compare Policy function F. u = -F*x
# Left term in F equal to 1 + left term in F_sol???
julia> F
1×2 Matrix{Float64}:
 0.866667  0.0

julia> F_sol
1×2 Matrix{Float64}:
 -0.133333  0.0

[andreasvarga]

The trick with (R+R')/2 works just in the case you indicated, but it does not work when solving the Riccati equation.

julia> using LinearAlgebra, MatrixEquations

julia> A = rand(2,2); B = rand(2,2); Q1 = rand(2,2); R1 = rand(2,2);

julia> Q = (Q1+Q1')/2;  R = (R1+R1')/2;

julia> X, EIGS, F = ared(A,B,R,Q);

julia> A'*X*A-X-A'*X*B*inv(R+B'*X*B)*B'*X*A+Q
2×2 Matrix{Float64}:
 -1.33227e-15  -2.77556e-16
  0.0           5.55112e-17

julia> A'*X*A-X-A'*X*B*inv(R1+B'*X*B)*B'*X*A+Q1
2×2 Matrix{Float64}:
 -0.0024168   0.147394
 -0.148686   -0.000172719

[baggepinnen]

But drawing random matrices is not guaranteed to yield positive definite matrices, the trick is only meant to make them symmetric, not positive definite.

[azev77]

I think I basically got it. But I had to derive things.
Bottom line: I need to replace S by 0.5*S b/c the formulation in this package is different from mine (in some way).
My Riccati equation has 0.5*S, your equation has S

using MatrixEquations, LinearAlgebra; 
δ=0.1; r=0.05; z= r + δ +.02; β=1.0/(1.0 + r); # parameters

I_SS = (z-r-δ)/(r+δ);
K_SS = I_SS/δ;
D_SS = z*K_SS - I_SS - 0.5*(I_SS)^(2.0);
V_SS = ((1.0+r)/r)*D_SS;
Q_SS = (1+r)*(1+I_SS);

P_sol = [(V_SS - Q_SS*K_SS) 0.0; Q_SS   0.0]; # closed form P. Value = x'*P*x
F_sol = [-I_SS   0.0];                        # closed form u. Policy = -F_sol*x

A=[1.0 0.0; 0.0  (1.0 - δ)];
B=[0.0; 1.0];
Q=[0.0 0.0; z 0.0];
R=[-0.5];
S=[-1.0; 0.0];

S = 0.5*S; # my Riccatti is different from MatrixEquations.jl
#Q = -1.0*Q; R = -1.0*R; S = -1.0*S;    # max to min. DO NOT NEED!
A = sqrt(β)*A; B = sqrt(β)*B; # scale by discount factor
Q = (Q + Q') ./ 2.0; R = (R + R') ./ 2.0;  P_sol = (P_sol + P_sol') ./ 2.0; # make symmetric  
                       
P, CLSEIG, F = ared(A,B,R,Q,S);             # solve for numerical value function = x'*P*x

# Compare Value function P. CORRECT!
julia> P
2×2 Matrix{Float64}:
 0.186667  0.595
 0.595     0.0

julia> P_sol
2×2 Matrix{Float64}:
 0.186667  0.595
 0.595     0.0

# Compare Policy function F. u = -F*x. CORRECT!
julia> F
1×2 Matrix{Float64}:
 -0.133333  -0.0

julia> F_sol
1×2 Matrix{Float64}:
 -0.133333  0.0

Fn  = pinv(R .+ (B')*P*B) * ((B')*P*A +S'); # numerical policy = -F_n*x. SAME as F from ared().
julia> Fn
1×2 Matrix{Float64}:
 -0.133333  0.0

PS: @baggepinnen was right
if I make the objective negative (since I'm maximizing)
Q = -1.0*Q; R = -1.0*R; S = -1.0*S;
The corresponding guess for the value function will also be negative -x' * P *x
So I will need P_sol = -1.0* P_sol
However, the minimizer F will be unchanged...

[azev77]

@andreasvarga
In continuous time I found I need to replace A with A = A - (ρ/2)*I and I get the correct answer.
Using A = A + (ρ/2)*I gives arec() an error...
My guess, it's prob a notational issue w/ how I set up the rest of my model...