@Juan, you ask about the future? I have no idea about the future, but I will show you two examples from the present where Julia easily beats C and Fortran (by C and Fortran I mean latest GCC compilers because languages themselves are irrelevant here). My dear friend, Julia toppled them all from the ring, still only one old stubborn player who refuses to die, Intel Fortran. The first example is a real-world one from A Comparison of Programming Languages in Economics and the second one is a toy example. I have many of these, but for the sake of brevity, I will show only two.
In the following examples, I will show Fortran and Julia codes. For compiling the Fortran code, all possible aggressive optimizations (for the GCC 8.1.0 compiler) were turned on via -Ofast and Julia is compiled with the -O3 --check-bounds=no options. Julia runs 2X faster and 1.5X faster in the examples below, respectively.
Here are the results first, and I will leave you with the code.
First example:
Elapsed time is 1.62500000 # Fortran
0.802962 seconds (1.11 k allocations: 353.815 MiB, 3.81% gc time) # Julia
Second example:
Time (sec): 0.5160000 # Fortran
0.352425 seconds (4 allocations: 160 bytes) # Julia
Now Fortran code (Example 1):
!============================================================================
! Name : RBC_F90.f90
! Description : Basic RBC model with full depreciation, Intel compiler
! Date : July 21, 2013
!============================================================================
program RBC_F90
!----------------------------------------------------------------
! 0. variables to be defined
!----------------------------------------------------------------
implicit none
integer, parameter :: nGridCapital = 17820
integer, parameter :: nGridProductivity = 5
real(8), parameter :: tolerance = 0.0000001
integer :: nCapital, nCapitalNextPeriod, gridCapitalNextPeriod, nProductivity, nProductivityNextPeriod
integer :: iteration
real :: start_time, stop_time
real(8) :: aalpha, bbeta, capitalSteadyState, outputSteadyState, consumptionSteadyState
real(8) :: valueHighSoFar, valueProvisional, consumption, capitalChoice
real(8) :: maxDifference,diff,diffHighSoFar
real(8), dimension(nGridProductivity) :: vProductivity
real(8), dimension(nGridProductivity,nGridProductivity) :: mTransition
real(8), dimension(nGridCapital) :: vGridCapital
real(8), dimension(nGridCapital,nGridProductivity) :: mOutput, mValueFunction, mValueFunctionNew, mPolicyFunction
real(8), dimension(nGridCapital,nGridProductivity) :: expectedValueFunction
!----------------------------------------------------------------
! 1. Calibration
!----------------------------------------------------------------
call cpu_time(start_time)
aalpha = 1.d0/3 ! Elasticity of output w.r.t
bbeta = 0.95d0 ! Discount factor
! Productivity value
vProductivity = (/0.9792, 0.9896, 1.0000, 1.0106, 1.0212/)
! Transition matrix
mTransition = reshape( (/0.9727, 0.0273, 0.0, 0.0, 0.0, &
0.0041, 0.9806, 0.0153, 0.0, 0.0, &
0.0, 0.0082, 0.9837, 0.0082, 0.0, &
0.0, 0.0, 0.0153, 0.9806, 0.0041, &
0.0, 0.0, 0.0, 0.0273, 0.9727 /), (/5,5/))
mTransition = transpose(mTransition)
!----------------------------------------------------------------
! 2. Steady state
!----------------------------------------------------------------
capitalSteadyState = (aalpha*bbeta)**(1.0/(1.0-aalpha))
outputSteadyState = capitalSteadyState**(aalpha)
consumptionSteadyState = outputSteadyState-capitalSteadyState
print *, 'Steady State values'
print *, 'Output: ', outputSteadyState, 'Capital: ', capitalSteadyState, 'Consumption: ', consumptionSteadyState
! Grid for capital
do nCapital = 1, nGridCapital
vGridCapital(nCapital) = 0.5*capitalSteadyState + 0.00001d0*(nCapital-1)
end do
!----------------------------------------------------------------
! 3. Pre-build Output for each point in the grid
!----------------------------------------------------------------
do nProductivity = 1, nGridProductivity
do nCapital = 1, nGridCapital
mOutput(nCapital, nProductivity) = vProductivity(nProductivity)*(vGridCapital(nCapital)**aalpha)
end do
end do
!----------------------------------------------------------------
! 4. Main Iteration
!----------------------------------------------------------------
maxDifference = 10.0
iteration = 0
do while (maxDifference > tolerance)
expectedValueFunction = matmul(mValueFunction,transpose(mTransition));
do nProductivity = 1,nGridProductivity
! We start from previous choice (monotonicity of policy function)
gridCapitalNextPeriod = 1
do nCapital = 1,nGridCapital
valueHighSoFar = -100000.0
do nCapitalNextPeriod = gridCapitalNextPeriod,nGridCapital
consumption = mOutput(nCapital,nProductivity)-vGridCapital(nCapitalNextPeriod)
valueProvisional = (1.0-bbeta)*log(consumption)+bbeta*expectedValueFunction(nCapitalNextPeriod,nProductivity)
if (valueProvisional>valueHighSoFar) then ! we break when we have achieved the max
valueHighSoFar = valueProvisional
capitalChoice = vGridCapital(nCapitalNextPeriod)
gridCapitalNextPeriod = nCapitalNextPeriod
else
exit
end if
end do
mValueFunctionNew(nCapital,nProductivity) = valueHighSoFar
mPolicyFunction(nCapital,nProductivity) = capitalChoice
end do
end do
maxDifference = maxval((abs(mValueFunctionNew-mValueFunction)))
mValueFunction = mValueFunctionNew
iteration = iteration + 1
! if (mod(iteration,10)==0 .OR. iteration==1) then
! print *, 'Iteration:', iteration, 'Sup Diff:', MaxDifference
! end if
end do
!----------------------------------------------------------------
! 5. PRINT RESULTS
!----------------------------------------------------------------
print *, 'Iteration:', iteration, 'Sup Diff:', MaxDifference, char(10)
print *, 'My check:', mPolicyFunction(1000,3)
call cpu_time(stop_time)
print *, 'Elapsed time is', stop_time - start_time
end program RBC_F90
! ==================================================== !
! ---------------------- OUTPUT ---------------------- !
! ==================================================== !
Steady State values
Output: 0.56273143387113778 Capital: 0.17819828739252699 Consumption: 0.38453314647861081
Iteration: 257 Sup Diff: 9.7160058332157462E-008
My check: 0.14654914369626351
Elapsed time is 1.62500000
Julia code (Example 1):
function main()
## 1. Calibration
α = 1/3 # Elasticity of output w.r.t. capital
β = 0.95 # Discount factor num_triangles
# 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]
mTransition = mTransition' # CHANGE: do this here, not inside the loop
# 2. Steady State
capitalSteadyState = (α*β)^(1.0/(1.0-α))
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 += 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
# =============================================== #
# -------------------- OUTPUT ------------------- #
# =============================================== #
julia> main()
julia> @time main()
Iteration = 257 Sup Diff = 9.716035653806188e-8
My check = 0.1465491436962635
My check = -0.9714880021887345
Iteration = 257 Sup Diff = 9.716035653806188e-8
My check = 0.1465491436962635
My check = -0.9714880021887345
0.802962 seconds (1.11 k allocations: 353.815 MiB, 3.81% gc time)
Fortran code (Example 2):
program test_prod
implicit none
integer, parameter :: n = 10**8
double precision, allocatable :: x(:), y(:)
integer :: t0, count_rate, count_max, t1
allocate ( x(n), y(n) )
call random_number(x)
call system_clock(t0,count_rate,count_max)
call toarr(x,y)
call system_clock(t1)
write (*,'(a,f9.7,a,5f9.5)') 'Time (sec): ',real(t1-t0)/count_rate, char(10), y(1:5)
contains
subroutine toarr(x,y)
double precision, intent(in) :: x(:)
double precision, intent(out):: y(:)
y = x*x
where (x <= 0.3) y = sqrt(x)
end subroutine
end program test_prod
! ==================================================== !
! ---------------------- OUTPUT ---------------------- !
! ==================================================== !
Time (sec): 0.5160000
0.19949 0.94330 0.11662 0.37159 0.15961
Julia code (Example 2):
function toarr!(x,y)
v = x[1]
for (i,v) in enumerate(x)
y[i] = ifelse(v > 0.3, v*v, sqrt(v))
end
end
# =============================================== #
# -------------------- OUTPUT ------------------- #
# =============================================== #
julia> x = rand(10^8);
julia> y = zeros(10^8);
julia> @time toarr!(x,y)
0.352425 seconds (4 allocations: 160 bytes)
julia> println(y[1:5])
[0.193422, 0.176242, 0.520951, 0.63757, 0.658653]