Hey all, I found this benchmark comparing different languages in economics and I wonder whether the underlying Julia code can be optimized further (full repository is available here):

```
## Basic RBC model with full depreciation
#
# Jesus Fernandez-Villaverde
# Haverford, July 29, 2013
function main()
## 1. Calibration
α = 1/3 # Elasticity of output w.r.t. capital
β = 0.95 # Discount factor
# Productivity values
vProductivity = [0.9792 0.9896 1.0000 1.0106 1.0212]
# Transition matrix
mTransition = [0.9727 0.0273 0.0000 0.0000 0.0000;
0.0041 0.9806 0.0153 0.0000 0.0000;
0.0000 0.0082 0.9837 0.0082 0.0000;
0.0000 0.0000 0.0153 0.9806 0.0041;
0.0000 0.0000 0.0000 0.0273 0.9727]
# 2. Steady State
capitalSteadyState = (α*β)^(1/(1-α))
outputSteadyState = capitalSteadyState^α
consumptionSteadyState = outputSteadyState-capitalSteadyState
println("Output = ",outputSteadyState," Capital = ",capitalSteadyState," Consumption = ",consumptionSteadyState)
# We generate the grid of capital
vGridCapital = collect(0.5*capitalSteadyState:0.00001:1.5*capitalSteadyState)
nGridCapital = length(vGridCapital)
nGridProductivity = length(vProductivity)
# 3. Required matrices and vectors
mOutput = zeros(nGridCapital,nGridProductivity)
mValueFunction = zeros(nGridCapital,nGridProductivity)
mValueFunctionNew = zeros(nGridCapital,nGridProductivity)
mPolicyFunction = zeros(nGridCapital,nGridProductivity)
expectedValueFunction = zeros(nGridCapital,nGridProductivity)
# 4. We pre-build output for each point in the grid
mOutput = (vGridCapital.^α)*vProductivity;
# 5. Main iteration
maxDifference = 10.0
tolerance = 0.0000001
iteration = 0
while(maxDifference > tolerance)
expectedValueFunction = mValueFunction*mTransition';
for nProductivity in 1:nGridProductivity
# We start from previous choice (monotonicity of policy function)
gridCapitalNextPeriod = 1
for nCapital in 1:nGridCapital
valueHighSoFar = -1000.0
capitalChoice = vGridCapital[1]
for nCapitalNextPeriod in gridCapitalNextPeriod:nGridCapital
consumption = mOutput[nCapital,nProductivity]-vGridCapital[nCapitalNextPeriod]
valueProvisional = (1-β)*log(consumption)+β*expectedValueFunction[nCapitalNextPeriod,nProductivity]
if (valueProvisional>valueHighSoFar)
valueHighSoFar = valueProvisional
capitalChoice = vGridCapital[nCapitalNextPeriod]
gridCapitalNextPeriod = nCapitalNextPeriod
else
break # We break when we have achieved the max
end
end
mValueFunctionNew[nCapital,nProductivity] = valueHighSoFar
mPolicyFunction[nCapital,nProductivity] = capitalChoice
end
end
maxDifference = maximum(abs.(mValueFunctionNew-mValueFunction))
mValueFunction, mValueFunctionNew = mValueFunctionNew, mValueFunction
iteration = iteration+1
if mod(iteration,10)==0 || iteration == 1
println(" Iteration = ", iteration, " Sup Diff = ", maxDifference)
end
end
println(" Iteration = ", iteration, " Sup Diff = ", maxDifference)
println(" ")
println(" My check = ", mPolicyFunction[1000,3])
println(" My check = ", mValueFunction[1000,3])
end
```

Right now, Julia 1.0 is around 40%-50% slower than both Fortran and C++ (in the paper, as well as on my machine).

As far as I can see, the code is written properly in column-major-order and I found no obvious type-instabilities. I also tried **Traceur** and `@code_warntype`

, but that didn’t get me more insights.

I have the highest optimization level enabled and removed some printing, so that Juno’s slow printing doesn’t mess up the results. I know that one could play around with some libraries for the matrix multiplication, but that will affect the results only marginally.

I would appreciate any tips!