Solving Bellman Equations by the Collocation Method

A large class of economic models involves solving for functional equations of the form:

equation1

A well known example is the stochastic optimal growth model. An agent owns a consumption good $y$ at time $t$, which can be consumed or invested. Next period output depends on how much is invested at time $t$ and on a shock $z$ realized at the end of the current period. One can think of a farmer deciding the quantity of seeds to be planted during the spring, taking into account weather forecast for the growing season.

A common technique for solving this class of problem is value function iteration. While value function iteration is quite intuitive (it follows the proof of existence of a solution to the functional equation above), it is not the only one available. This post describes the collocation method, which transforms the problem of finding a function into a root-finding one. The gain from this change of perspective is that any root-finding algorithm can then be used. In particular, one may use the Newton method, which converges at a quadratic rate in the neighborhood of the solution if the function is smooth enough.

Value function iteration

Value function iteration takes advantage of the fact that the Bellman operator $T$ is a contraction mapping on the set of continuous bounded functions on $\mathbb R_+$ under the supremum distance

equation2

An immediate consequence if that the sequence $w,Tw,T^2w$,… converges uniformly to $w$ (starting with with any bounded and continuous $w$). The following code in Julia v0.6 illustrates the convergence of the series ${T^nw}$.

 #=
Julien Pascal

Code heavily based on:
----------------------
https://lectures.quantecon.org/jl/optgrowth.html
by Spencer Lyon, John Stachurski

I have only made minor modifications:
#------------------------------------
* I added a type optGrowth
* I use the package Interpolations
* I calculate the expectation w.r.t the aggregate
  shock using a Gauss-Legendre quadrature scheme
  instead of Monte-Carlo


=#

using QuantEcon
using Optim
using CompEcon
using PyPlot
using Interpolations
using FileIO


type optGrowth

  w::Array{Float64,1}
  β::AbstractFloat
  grid::Array{Float64,1}
  u::Function
  f::Function
  shocks::Array{Float64,1}
  Tw::Array{Float64,1}
  σ::Array{Float64,1}
  el_k::Array{Float64,1}
  wl_k::Array{Float64,1}
  compute_policy::Bool
  w_func::Function

end

function optGrowth(;w = Array{Float64,1}[],
                    α = 0.4,
                    β = 0.96,
                    μ = 0,
                    s = 0.1,
                    grid_max = 4,         # Largest grid point
                    grid_size = 200,      # Number of grid points
                    shock_size = 250,     # Number of shock draws in Monte Carlo integral
                    Tw = Array{Float64,1}[],
                    σ = Array{Float64,1}[],
                    el_k = Array{Float64,1}[],
                    wl_k = Array{Float64,1}[],
                    compute_policy = true
                  )


  grid_y = collect(linspace(1e-5, grid_max, grid_size))
  shocks = exp.(μ + s * randn(shock_size))

  # Utility
  u(c) = log(c)
  # Production
  f(k) = k^α

  w = 5 * log.(grid_y)
  Tw = 5 * log.(grid_y)
  σ = 5 * log.(grid_y)

  el_k, wl_k = qnwlogn(10, μ, s^2) #10 weights and nodes for LOG(e_t) distributed N(μ,s^2)

  w_func = x -> x
  optGrowth(
    w,
    β,
    grid_y,
    u,
    f,
    shocks,
    Tw,
    σ,
    el_k,
    wl_k,
    compute_policy,
    w_func
    )
end

"""
The approximate Bellman operator, which computes and returns the
updated value function Tw on the grid points.

#### Arguments

`model` : a model of type optGrowth

`Modifies model.σ, model.w and model.Tw

"""
function bellman_operator!(model::optGrowth)

    # === Apply linear interpolation to w === #
    knots = (model.grid,)
    itp = interpolate(knots, model.w, Gridded(Linear()))

    #w_func(x) = itp[x]

    model.w_func = x -> itp[x]

    if model.compute_policy
        model.σ = similar(model.w)
    end

    # == set Tw[i] = max_c { u(c) + β E w(f(y  - c) z)} == #
    for (i, y) in enumerate(model.grid)

        #Monte Carlo
        #-----------
        #objective(c) = - model.u(c) - model.β * mean(w_func(model.f(y - c) .* model.shocks))

        #Gauss-Legendre
        #--------------
        function objective(c)

          expectation = 0.0

          for k = 1:length(model.wl_k)
            expectation += model.wl_k[k]*(model.w_func(model.f(y - c) * model.el_k[k]))
          end

          - model.u(c) - model.β * expectation

        end

        res = optimize(objective, 1e-10, y)

        if model.compute_policy
            model.σ[i] = Optim.minimizer(res)
        end

        model.Tw[i] = - Optim.minimum(res)
        model.w[i] = - Optim.minimum(res)
    end


end


model = optGrowth()

function solve_optgrowth!(model::optGrowth;
                         tol::AbstractFloat=1e-6,
                         max_iter::Integer=500)

    w_old = copy(model.w)  # Set initial condition
    error = tol + 1
    i = 0


    # Iterate to find solution
    while i < max_iter

        #update model.w
        bellman_operator!(model)

        error = maximum(abs, model.w - w_old)

        if error < tol
          break
        end

        w_old = copy(model.w)
        i += 1
    end

end

#-----------------------------------
# Solve by value function iteration
#-----------------------------------
@time solve_optgrowth!(model)
# 3.230501 seconds (118.18 M allocations: 1.776 GiB, 3.51% gc time)


#-------------------------------
# Compare with the true solution
#-------------------------------
α = 0.4
β = 0.96
μ = 0
s = 0.1

c1 = log(1 - α * β) / (1 - β)
c2 = (μ + α * log(α * β)) / (1 - α)
c3 = 1 / (1 - β)
c4 = 1 / (1 - α * β)


# True optimal policy
c_star(y) = (1 - α * β) * y

# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log.(y)

fig, ax = subplots(figsize=(9, 5))
ax[:set_ylim](-35, -24)
ax[:plot](model.grid, model.w_func.(model.grid), lw=2, alpha=0.6, label="approximate value function")
ax[:plot](model.grid, v_star.(model.grid), lw=2, alpha=0.6, label="true value function")
ax[:legend](loc="lower right")

fig, ax = subplots(figsize=(9, 5))
ax[:set_xlim](0.1, 4.0)
ax[:set_ylim](0.00, 0.2)
ax[:plot](model.grid, abs.(model.w_func.(model.grid) -v_star.(model.grid)), lw=2, alpha=0.6, label="error")
ax[:legend](loc="lower right")



Output

VFI VFI

The collocation method

The collocation method takes a different route. Let us remember that we are looking for a function $w$. Instead of solving for the values of $w$ on a grid and then interpolating, why not looking for a function directly? To do so, let us assume that $w$ can reasonably be approximated by a function $\hat{w}$:

equation3

with $ \phi_1(x) $ , $ \phi_2(x) $,…, $ \phi_n(x) $ a set of linearly independent basis functions and $c_1$, $c_2$, …, $c_n$ $n$ coefficient to be found. Replacing $w(x)$ with $\hat{w(x)}$ into the functional equation and reorganizing gives:

equation4

This equation has to hold (almost) exactly at $n$ points (also called nodes): $y_1$, $y_2$, …, $y_n$:

equation5

The equation above defines a system of $n$ equation with as many unknown, which can be compactly written as: $$ f(\boldsymbol{c}) = \boldsymbol{0} $$

Newton or quasi-Newton can be used to solve for the root of $f$. In the code that follows, I use Broyden’s method. Let us illustrate this technique using a Chebychev polynomial basis and Chebychev nodes. In doing so, we avoid Runge’s phenomenon associated with a uniform grid.

Implementation using CompEcon

 #---------------------------------------------

 # Julien Pascal

 # Solve the stochastic optimal growth problem
 # using the collocation method

 #---------------------------------------------




 using QuantEcon
 using Optim
 using CompEcon
 using PyPlot
 using Interpolations



 type optGrowthCollocation

   w::Array{Float64,1}
   β::AbstractFloat
   grid::Array{Float64,1}
   u::Function
   f::Function
   shocks::Array{Float64,1}
   Tw::Array{Float64,1}
   σ::Array{Float64,1}
   el_k::Array{Float64,1}
   wl_k::Array{Float64,1}
   compute_policy::Bool
   order_approximation::Int64 #number of element in the functional basis along each dimension
   functional_basis_type::String #type of functional basis
   fspace::Dict{Symbol,Any} #functional basis
   fnodes::Array{Float64,1} #collocation nodes
   residual::Array{Float64,1} #vector of residual. Should be close to zero
   a::Array{Float64,1} #polynomial coefficients
   w_func::Function

 end


   #####################################
   # Function that finds a solution
   # to f(x) = 0
   # using Broyden's "good" method
   # and a simple backstepping procedure as described
   # in Miranda and Fackler (2009)
   #
   # input :
   # --------
   # * x0:                 initial guess for the root
   # * f:                  function in f(x) = 0
   # * maxit:              maximum number of iterations
   # * tol:                tolerance level for the zero
   # * fjavinc:            initial inverse of the jacobian. If not provided, then inverse of the
   #                       Jacobian is calculated by finite differences
   # * maxsteps:           maximum number of backsteps
   # * recaculateJacobian: number of iterations in-between two calculations of the Jacobian
   #
   # output :
   # --------
   # * x: one zero of f
   # * it: number of iterations necessary to reached the solution
   # * fjacinv: pseudo jacobian at the last iteration
   # * fnorm: norm f(x) at the last iteration
   #
   #######################################
   function find_broyden(x0::Vector, f::Function, maxit::Int64, tol::Float64, fjacinv = eye(length(x0));
                         maxsteps = 5, recaculateJacobian = 1)

       println("a0 = $(x0)")
       fnorm = tol*2
       it2 = 0 #to re-initialize the jacobian

       ################################
       #initialize guess for the matrix
       ################################
       fjacinv_function = x-> Calculus.finite_difference_jacobian(f, x)
       #fjacinv_function = x -> ForwardDiff.gradient(f, x)

       # If the user do not provide an initial guess for the jacobian
       # One is calculated using finite differences.
       if fjacinv == eye(length(x0))
           ################################################
           # finite differences to approximate the Jacobian
           # at the initial value
           # this is slow. Seems to improve performances
           # when x0 is of high dimension.
           println("Calculating the Jacobian by finite differences")
           #@time fjacinv = Calculus.finite_difference_jacobian(f, x0)
           @time fjacinv = fjacinv_function(x0)

           println("Inverting the Jacobian")
           try
               fjacinv = inv(fjacinv)
           catch
               try
                   println("Jacobian non-invertible\n calculating pseudo-inverse")
                   fjacinv = pinv(A)
               catch
                   println("Failing Calculating the pseudo-inverse. Initializing with In")
                   fjacinv = eye(length(x0))
               end
           end
           println("Done")
       else
           println("Using User's input as a guess for the Jacobian.")
       end

       fval = f(x0)

       for it=1:maxit

           it2 +=1

           #every 30 iterations, reinitilize the jacobian
           if mod(it2, recaculateJacobian) == 0

               println("Re-calculating the Jacobian")

               fjacinv = fjacinv_function(x0)

               try
                   fjacinv = inv(fjacinv)
               catch
                   try
                       println("Jacobian non-invertible\n calculating pseudo-inverse")
                       fjacinv = pinv(A)
                   catch
                       println("Failing Calculating the pseudo-inverse. Initializing with In")
                       fjacinv = eye(length(x0))
                   end
               end

           end


           println("it = $(it)")

           fnorm = norm(fval)

           if fnorm < tol
               println("fnorm = $(fnorm)")
               return x0, it, fjacinv, fnorm
           end

           d = -(fjacinv*fval)

           fnormold = Inf
           ########################
           # Backstepping procedure
           ########################
           for backstep = 1:maxsteps

               if backstep > 1
                   println("backstep = $(backstep-1)")
               end

               fvalnew = f(x0 + d)
               fnormnew = norm(fvalnew)

               if fnormnew < fnorm
                   break
               end

               if fnormold < fnormnew
                   d=2*d
                   break
               end

               fnormold = fnormnew

               d = d/2

           end
           ####################
           ####################

           x0 = x0 + d

           fold = fval
           fval = f(x0)

           u = fjacinv*(fval - fold)

           #Update the pseudo Jacobian:
           fjacinv = fjacinv + ((d-u)*(transpose(d)*fjacinv))/(dot(d,u))

           println("a$(it) = $(x0)")
           println("fnorm = $(fnorm)")

           if isnan.(x0) == trues(length(x0))
               println("Error. a$(it) = NaN for each component")
               x0 = zeros(length(x0))
               return x0, it, fjacinv, fnorm
           end
       end

       println("In function find_broyden\n")
       println("Maximum number of iterations reached.\n")
       println("No convergence.")
       println("Returning fnorm = NaN as a solution")
       fnorm = NaN
       return x0, maxit, fjacinv, fnorm

   end

 function optGrowthCollocation(;w = Array{Float64,1}[],
                               α = 0.4,
                               β = 0.96,
                               μ = 0,
                               s = 0.1,
                               grid_max = 4,         # Largest grid point
                               grid_size = 200,      # Number of grid points
                               shock_size = 250,     # Number of shock draws in Monte Carlo integral
                               Tw = Array{Float64,1}[],
                               σ = Array{Float64,1}[],
                               el_k = Array{Float64,1}[],
                               wl_k = Array{Float64,1}[],
                               compute_policy = true,
                               order_approximation = 40,
                               functional_basis_type = "chebychev",
                             )


   grid_y = collect(linspace(1e-5, grid_max, grid_size))
   shocks = exp.(μ + s * randn(shock_size))

   # Utility
   u(c) = log.(c)
   # Production
   f(k) = k^α

   el_k, wl_k = qnwlogn(10, μ, s^2) #10 weights and nodes for LOG(e_t) distributed N(μ,s^2)

   lower_bound_support = minimum(grid_y)
   upper_bound_support = maximum(grid_y)

   n_functional_basis = [order_approximation]

   if functional_basis_type == "chebychev"
       fspace = fundefn(:cheb, n_functional_basis, lower_bound_support, upper_bound_support)
   elseif functional_basis_type == "splines"
       fspace = fundefn(:spli, n_functional_basis, lower_bound_support, upper_bound_support, 1)
   elseif functional_basis_type == "linear"
       fspace = fundefn(:lin, n_functional_basis, lower_bound_support, upper_bound_support)
   else
       error("functional_basis_type has to be either chebychev, splines or linear.")
   end


   fnodes = funnode(fspace)[1]
   residual = zeros(size(fnodes)[1])
   a = ones(size(fnodes)[1])

   w = ones(size(fnodes)[1])
   Tw = ones(size(fnodes)[1])
   σ = ones(size(fnodes)[1])

   w_func = x-> x

   optGrowthCollocation(
     w,
     β,
     grid_y,
     u,
     f,
     shocks,
     Tw,
     σ,
     el_k,
     wl_k,
     compute_policy,
     order_approximation,
     functional_basis_type,
     fspace,
     fnodes,
     residual,
     a,
     w_func
     )
 end



 function residual!(model::optGrowthCollocation)


     model.w_func = y -> funeval(model.a, model.fspace, [y])[1][1]   


     function objective(c, y)

       expectation = 0.0

       for k = 1:length(model.wl_k)
         expectation += model.wl_k[k]*(model.w_func(model.f(y - c) * model.el_k[k]))
       end

       - model.u(c) - model.β * expectation

     end

     # Loop over nodes
     for i in 1:size(model.fnodes)[1]

         y = model.fnodes[i,1]

         res = optimize(c -> objective(c, y), 1e-10, y)

         if model.compute_policy
             model.σ[i] = Optim.minimizer(res)
         end

         model.Tw[i] = - Optim.minimum(res)
         model.w[i] = model.w_func(y)

         model.residual[i] = - model.w[i] + model.Tw[i]
     end

 end


 model = optGrowthCollocation(functional_basis_type = "chebychev")

 residual!(model)

 function solve_optgrowth!(model::optGrowthCollocation;
                          tol=1e-6,
                          max_iter=500)

     # Initialize guess for coefficients
     # by giving the "right shape"
     # ---------------------------------
     function objective_initialize!(x, model)

       #update polynomial coeffficients
       model.a = copy(x)

       model.w_func = y -> funeval(model.a, model.fspace, [y])[1][1]

       return abs.(model.w_func.(model.fnodes[:,1]) - 5.0 * log.(model.fnodes[:,1]))

     end


     minx, iterations, Jac0, fnorm = find_broyden(model.a, x -> objective_initialize!(x, model), max_iter, tol)


     # Solving the model by collocation
     # using the initial guess calculated above
     #-----------------------------------------
     function objective_residual!(x, model)

       #update polynomial coeffficients
       model.a = copy(x)

       #calculate residual
       residual!(model)

       return abs.(model.residual)

     end

     minx, iterations, Jac, fnorm = find_broyden(model.a, x -> objective_residual!(x, model), max_iter, tol)


 end

 #-----------------------------------
 # Solve by collocation
 #-----------------------------------
 @time solve_optgrowth!(model)


 #-------------------------------
 # Compare with the true solution
 #-------------------------------
 α = 0.4
 β = 0.96
 μ = 0
 s = 0.1

 c1 = log(1 - α * β) / (1 - β)
 c2 = (μ + α * log(α * β)) / (1 - α)
 c3 = 1 / (1 - β)
 c4 = 1 / (1 - α * β)


 # True optimal policy
 c_star(y) = (1 - α * β) * y

 # True value function
 v_star(y) = c1 + c2 * (c3 - c4) + c4 * log.(y)

 fig, ax = subplots(figsize=(9, 5))
 ax[:set_ylim](-35, -24)
 ax[:plot](model.grid, model.w_func.(model.grid), lw=2, alpha=0.6, label="approximate value function")
 ax[:plot](model.grid, v_star.(model.grid), lw=2, alpha=0.6, label="true value function")
 ax[:legend](loc="lower right")


 fig, ax = subplots(figsize=(9, 5))
 ax[:set_xlim](0.1, 4.0)
 ax[:set_ylim](0.01, 0.2)
 ax[:plot](model.grid, abs.(model.w_func.(model.grid) -v_star.(model.grid)), lw=2, alpha=0.6, label="error")
 ax[:legend](loc="lower right")


Output

VFI VFI

We obtain more are less the same function $w$. But it turns out that the value function iteration implementation is much faster. One reason seems to be the efficiency associated with the package Interpolations: it is more than 20 times faster to evaluate $w$ using the package Interpolations rather than using the package CompEcon:

  #using Interpolations
  #--------------------
  @time for i=1:1000000
   model.w_func.(model.grid[1])
  end
  #0.230861 seconds (2.00 M allocations: 30.518 MiB, 1.28% gc time)

  #using CompEcon
  #--------------
  @time for i=1:1000000
   model.w_func.(model.grid[1])
  end
  # 4.998902 seconds (51.00 M allocations: 3.546 GiB, 13.39% gc time)

Implementation using ApproxFun

Significant speed-up can be obtained by using the package ApproxFun, as illustrated by the code below.


#---------------------------------------------

# Julien Pascal

# Solve the stochastic optimal growth problem
# using the collocation method

# Implementation using ApproxFun

#---------------------------------------------




using QuantEcon
using Optim
using CompEcon
using PyPlot
using Interpolations
using FileIO
using ApproxFun
using ProfileView



type optGrowthCollocation

  w::Array{Float64,1}
  β::AbstractFloat
  grid::Array{Float64,1}
  u::Function
  f::Function
  shocks::Array{Float64,1}
  Tw::Array{Float64,1}
  σ::Array{Float64,1}
  el_k::Array{Float64,1}
  wl_k::Array{Float64,1}
  compute_policy::Bool
  order_approximation::Int64 #number of element in the functional basis along each dimension
  functional_basis_type::String #type of functional basis
  fspace::Dict{Symbol,Any} #functional basis
  fnodes::Array{Float64,1} #collocation nodes
  residual::Array{Float64,1} #vector of residual. Should be close to zero
  a::Array{Float64,1} #polynomial coefficients
  fApprox::ApproxFun.Fun{ApproxFun.Chebyshev{ApproxFun.Segment{Float64},Float64},Float64,Array{Float64,1}}
  w_func::Function
  tol::Float64

end


  #####################################
  # Function that find a solution
  # to f(x) = 0
  # using Broyden's "good" method
  # and simple backstepping procedure as described
  # in Miranda and Fackler (2009)
  #
  # input :
  # --------
  # * x0:                 initial guess for the root
  # * f:                  function in f(x) = 0
  # * maxit:              maximum number of iterations
  # * tol:                tolerance level for the zero
  # * fjavinc:            initial inverse of the jacobian. If not provided, then inverse of the
  #                       Jacobian is calculated by finite differences
  # * maxsteps:           maximum number of backsteps
  # * recaculateJacobian: number of iterations in-between two calculations of the Jacobian
  #
  # output :
  # --------
  # * x: one zero of f
  # * it: number of iterations necessary to reached the solution
  # * fjacinv: pseudo jacobian at the last iteration
  # * fnorm: norm f(x) at the last iteration
  #
  #######################################
  function find_broyden(x0::Vector, f::Function, maxit::Int64, tol::Float64, fjacinv = eye(length(x0));
                        maxsteps = 5, recaculateJacobian = 1)

      println("a0 = $(x0)")
      fnorm = tol*2
      it2 = 0 #to re-initialize the jacobian

      ################################
      #initialize guess for the matrix
      ################################
      # with Calculus
      #--------------
      fjacinv_function = x-> Calculus.finite_difference_jacobian(f, x)


      # If the user do not provide an initial guess for the jacobian
      # One is calculated using finite differences.
      if fjacinv == eye(length(x0))
          ################################################
          # finite differences to approximate the Jacobian
          # at the initial value
          # this is slow. Seems to improve performances
          # when x0 is of high dimension.
          println("Calculating the Jacobian by finite differences")
          #@time fjacinv = Calculus.finite_difference_jacobian(f, x0)
          @time fjacinv = fjacinv_function(x0)

          println("Inverting the Jacobian")
          try
              fjacinv = inv(fjacinv)
          catch
              try
                  println("Jacobian non-invertible\n calculating pseudo-inverse")
                  fjacinv = pinv(A)
              catch
                  println("Failing Calculating the pseudo-inverse. Initializing with In")
                  fjacinv = eye(length(x0))
              end
          end
          println("Done")
      else
          println("Using User's input as a guess for the Jacobian.")
      end

      fval = f(x0)

      for it=1:maxit

          it2 +=1

          #every 30 iterations, reinitilize the jacobian
          if mod(it2, recaculateJacobian) == 0

              println("Re-calculating the Jacobian")

              fjacinv = fjacinv_function(x0)

              try
                  fjacinv = inv(fjacinv)
              catch
                  try
                      println("Jacobian non-invertible\n calculating pseudo-inverse")
                      fjacinv = pinv(A)
                  catch
                      println("Failing Calculating the pseudo-inverse. Initializing with In")
                      fjacinv = eye(length(x0))
                  end
              end

          end


          println("it = $(it)")

          fnorm = norm(fval)

          if fnorm < tol
              println("fnorm = $(fnorm)")
              return x0, it, fjacinv, fnorm
          end

          d = -(fjacinv*fval)

          fnormold = Inf
          ########################
          # Backstepping procedure
          ########################
          for backstep = 1:maxsteps

              if backstep > 1
                  println("backstep = $(backstep-1)")
              end

              fvalnew = f(x0 + d)
              fnormnew = norm(fvalnew)

              if fnormnew < fnorm
                  break
              end

              if fnormold < fnormnew
                  d=2*d
                  break
              end

              fnormold = fnormnew

              d = d/2

          end
          ####################
          ####################

          x0 = x0 + d

          fold = fval
          fval = f(x0)

          u = fjacinv*(fval - fold)

          #Update the pseudo Jacobian:
          fjacinv = fjacinv + ((d-u)*(transpose(d)*fjacinv))/(dot(d,u))

          println("a$(it) = $(x0)")
          println("fnorm = $(fnorm)")

          if isnan.(x0) == trues(length(x0))
              println("Error. a$(it) = NaN for each component")
              x0 = zeros(length(x0))
              return x0, it, fjacinv, fnorm
          end
      end

      println("In function find_broyden\n")
      println("Maximum number of iterations reached.\n")
      println("No convergence.")
      println("Returning fnorm = NaN as a solution")
      fnorm = NaN
      return x0, maxit, fjacinv, fnorm

  end

function optGrowthCollocation(;w = Array{Float64,1}[],
                              α = 0.4,
                              β = 0.96,
                              μ = 0,
                              s = 0.1,
                              grid_max = 4,         # Largest grid point
                              grid_size = 200,      # Number of grid points
                              shock_size = 250,     # Number of shock draws in Monte Carlo integral
                              Tw = Array{Float64,1}[],
                              σ = Array{Float64,1}[],
                              el_k = Array{Float64,1}[],
                              wl_k = Array{Float64,1}[],
                              compute_policy = true,
                              order_approximation = 40,
                              functional_basis_type = "chebychev",
                            )


  grid_y = collect(linspace(1e-4, grid_max, grid_size))
  shocks = exp.(μ + s * randn(shock_size))

  # Utility
  u(c) = log.(c)
  # Production
  f(k) = k^α

  el_k, wl_k = qnwlogn(10, μ, s^2) #10 weights and nodes for LOG(e_t) distributed N(μ,s^2)

  lower_bound_support = minimum(grid_y)
  upper_bound_support = maximum(grid_y)

  n_functional_basis = [order_approximation]

  if functional_basis_type == "chebychev"
      fspace = fundefn(:cheb, n_functional_basis, lower_bound_support, upper_bound_support)
  else
      error("functional_basis_type has to be \"chebychev\" ")
  end


  fnodes = funnode(fspace)[1]
  residual = zeros(size(fnodes)[1])
  a = ones(size(fnodes)[1])

  tol = 0.001

  fApprox = (Fun(Chebyshev((minimum(grid_y))..(maximum(grid_y))), a))
  #fApprox = (Fun(Chebyshev(0..maximum(model.grid)), a))

  w_func = x-> fApprox(x)

  w = ones(size(fnodes)[1])
  Tw = ones(size(fnodes)[1])
  σ = ones(size(fnodes)[1])



  optGrowthCollocation(
    w,
    β,
    grid_y,
    u,
    f,
    shocks,
    Tw,
    σ,
    el_k,
    wl_k,
    compute_policy,
    order_approximation,
    functional_basis_type,
    fspace,
    fnodes,
    residual,
    a,
    fApprox,
    w_func,
    tol
    )
end



function residual!(model::optGrowthCollocation)


    model.fApprox = (Fun(Chebyshev((minimum(model.grid))..(maximum(model.grid))), model.a))
    model.w_func = x-> model.fApprox(x)

    function objective(c, y)

      expectation = 0.0

      for k = 1:length(model.wl_k)
        expectation += model.wl_k[k]*(model.w_func(model.f(y - c) * model.el_k[k]))
      end

      - model.u(c) - model.β * expectation

    end

    # Loop over nodes
    for i in 1:size(model.fnodes)[1]

        y = model.fnodes[i,1]

        res = optimize(c -> objective(c, y), 1e-10, y)

        if model.compute_policy
            model.σ[i] = Optim.minimizer(res)
        end

        model.Tw[i] = - Optim.minimum(res)
        model.w[i] = model.w_func(y)

        model.residual[i] = - model.w[i] + model.Tw[i]
    end

end


function solve_optgrowth!(model::optGrowthCollocation;
                         tol=1e-6,
                         max_iter=500)

    # Initialize guess for coefficients
    # by giving the "right shape"
    # ---------------------------------
    function objective_initialize!(x, model)

      #update polynomial coeffficients
      model.a = copy(x)

      model.fApprox = (Fun(Chebyshev((minimum(model.grid))..(maximum(model.grid))), model.a))

      model.w_func = x-> model.fApprox(x)

      return abs.(model.w_func.(model.fnodes[:,1]) - 5.0 * log.(model.fnodes[:,1]))

    end


    minx, iterations, Jac0, fnorm = find_broyden(model.a, x -> objective_initialize!(x, model), max_iter, tol)


    # Solving the model by collocation
    # using the initial guess calculated above
    #-----------------------------------------
    function objective_residual!(x, model)

      #update polynomial coeffficients
      model.a = copy(x)

      #calculate residual
      residual!(model)

      return abs.(model.residual)

    end

    minx, iterations, Jac, fnorm = find_broyden(model.a, x -> objective_residual!(x, model), max_iter, tol,  recaculateJacobian = 1)


end

#-----------------------------------
# Solve by collocation
#-----------------------------------
model = optGrowthCollocation(functional_basis_type = "chebychev")
@time solve_optgrowth!(model)
# 15.865923 seconds (329.12 M allocations: 4.977 GiB, 5.55% gc time)


#-------------------------------
# Compare with the true solution
#-------------------------------
α = 0.4
β = 0.96
μ = 0
s = 0.1

c1 = log(1 - α * β) / (1 - β)
c2 = (μ + α * log(α * β)) / (1 - α)
c3 = 1 / (1 - β)
c4 = 1 / (1 - α * β)


# True optimal policy
c_star(y) = (1 - α * β) * y

# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log.(y)

fig, ax = subplots(figsize=(9, 5))
ax[:set_ylim](-35, -24)
ax[:plot](model.grid, model.w_func.(model.grid), lw=2, alpha=0.6, label="approximate value function")
ax[:plot](model.grid, v_star.(model.grid), lw=2, alpha=0.6, label="true value function")

fig, ax = subplots(figsize=(9, 5))
ax[:set_xlim](0.1, 4.0)
ax[:set_ylim](- 0.05, 0.05)
ax[:plot](model.grid, abs.(model.w_func.(model.grid) -v_star.(model.grid)), lw=2, alpha=0.6, label="error")
ax[:legend](loc="lower right")


Output

VFI VFI

Computing time is divided by approximately $5$ using ApproxFun rather than CompEcon. But the value function iteration implementation is still the fastest one. One bottleneck seems to be the calculation of the Jacobian by finite differences. It is likely that using automatic differentiation would result is significant gain of time.

@time for i=1:1000000
 model.w_func.(model.grid[1])
end
# 0.391440 seconds (2.00 M allocations: 30.518 MiB, 1.37% gc time)