Hi! I am trying to push this code to be as efficient as possible. If I could achieve any performance gain it would be great because I am calling this function thousands of time.
The problem is essentially a problem of backward induction problem. In the last period, I need to choose the option that maximizes some payoff and in the periods before I need to maximize a payoff knowing that in the future I am behaving “optimally”.
If you consider there is a more efficient solution “paradigm” or something aside from small perfomance tweaks I would be super interested as well.
using BenchmarkTools
using ProgressMeter
using BitOperations
using Distributions
using LinearAlgebra
J = 1000
T = 200
Π = 100*rand(J,T)
Dᵢ = rand([0],J,1)
feat1 = rand([1,2,3,4,5,6,7,8],J,1)
feat2 = rand([1,2,3,4,5,6,7,8,9,10],J,1)
Φ1 = 10
Φ2 = 10
Sj = 20*rand(J,T)
Fj = 20*rand(J,T)
Π=Π.-Sj.-Fj
# This is the Backward Induction Algorithm with independent choices
@views function myfun2(Π,Fj,Sj,prob,SeqShocks,Dᵢ,β,J,T)
D = zeros(Int64,2,2,T,J)
v = zeros(Float64,2,2,T,J)
Dstar = zeros(Int64,J,T+1)
Vstar = zeros(Float64,J,T+1)
Dstar[:,1]=Dᵢ[1:J]
S = [0 1]
Shock=[0 999999999999999999999999999]
func2inner1!(Π, Sj, S, D, v, prob, Shock, β, J, T)
funct2inner2!(Dstar, Vstar, D, v, SeqShocks, T, J)
return Dstar[:,2:end], Vstar[:,2:end]
end
function func2inner1!(Π, Sj, S, D, v, prob, Shock, β, J, T)
Ev = zeros(Float64,2,T,J)
@inbounds for t ∈ 0:T-1, j ∈ 1:J, i ∈ 1:2
for ii ∈ 1:2
if t==0 && ii == 1
payoff = Π[j,T-t]+Sj[j,T-t]*S[i]-Shock[ii]
D[i,ii,T-t,j] = (payoff>0)
v[i,ii,T-t,j] = D[i,ii,T-t,j]*payoff
elseif ii==1
payoff_enter = Π[j,T-t]+Sj[j,T-t]*S[i]+β*Ev[2,T-t+1,j]
payoff_exit = β*Ev[1,T-t+1,j]
d=(payoff_enter>payoff_exit)
D[i,ii,T-t,j] = d
v[i,ii,T-t,j] = d*(payoff_enter)+(1-d)*payoff_exit
else
D[i,ii,T-t,j] = 0
payoff_exit = β*Ev[1,T-t+1,j]
v[i,ii,T-t,j] = payoff_exit
end
end
Ev[i,T-t,j]=prob*v[i,1,T-t,j]+(1-prob)*v[i,2,T-t,j]
end
end
function funct2inner2!(Dstar, Vstar, D, v, SeqShocks, T, J)
@inbounds for j=1:J, t=2:T+1
Dstar[j,t]=D[Dstar[j,t-1]+1,SeqShocks[j,t-1],t-1,j]
Vstar[j,t]=v[Dstar[j,t-1]+1,SeqShocks[j,t-1],t-1,j]
end
end
println("TIMES")
p = 0.9
SeqShocks = (rand(J,T).>p).+1
BI = @btime myfun2(Π,Fj,Sj,p,SeqShocks,Dᵢ,0.9,J,T)