jumpsolver.jl

using JuMP, Ipopt
m = Model(solver = IpoptSolver())

# Parameters
const γ = 2
const ρ = 0.05
const r = 0.041
const μR = 0.051
const ζ = 1.5
const σR = sqrt((ζ/γ+ 1)/2*(μR-r)^2/(ρ - r))
const n = 5000
const amin = -0.3
const amax = 1e4
x = map(x-> x + 5 * x^10, linspace(0, 1, n))
a = (amax-amin)/(maximum(x) - minimum(x))*x + amin
z = 3.0 * [0.01, 0.03]
λ = [0.5, 0.5]


# Variables to minimize
@variable(m, V[i=1:n, j=1:2], start = (z[j] +  r * a[i])^(1-γ) / (1-γ) / ρ)

# Expressions
## Constants
## First Derivative
### Backward
@NLexpression(m, Vpb1[j=1:2], (z[j] + a[1] * r)^(-γ))
@NLexpression(m, Vpbo[i=2:n, j=1:2], (V[i,j] - V[i-1,j]) / (a[i] - a[i-1]))
@NLexpression(m, Vpb[i=1:n, j=1:2], (i==1) * Vpb1[j] + (i>1) * Vpbo[max(i, 2), j])

### Forward
@NLexpression(m, Vpfo[i=1:(n-1), j=1:2], (V[i+1,j] - V[i,j]) / (a[i+1] - a[i]))
@NLexpression(m, Vpfn[j=1:2], (z[j] + a[n] * r + (μR-r)^2/(γ * σR^2) * a[n])^(-γ))
@NLexpression(m, Vpf[i=1:n, j=1:2], (i==n) * Vpfn[j] + (i<n) * Vpfo[min(i, n-1), j])

### All
@NLexpression(m, Vp[i=1:n, j=1:2, s=1:2], (s==1) * Vpb[i, j] + (s==2) * Vpf[i, j])

## Second Derivative
@NLexpression(m, Vpp1[j=1:2], 0.0)
@NLexpression(m, Vppn[j=1:2], - γ * Vp[n, j, 1] / a[n])
@NLexpression(m, Vppo[i=2:(n-1), j=1:2], ((a[i] - a[i-1]) * V[i+1,j] + (a[i+1] - a[i]) * V[i-1,j] - (a[i+1]- a[i-1]) * V[i,j]) / (0.5 * (a[i+1] - a[i-1]) * (a[i+1] - a[i]) * (a[i] - a[i-1])))
@NLexpression(m, Vpp[i=1:n, j=1:2], (i==1) * Vpp1[j] + (i==n) * Vppn[j] + (i>1) * (i<n) * Vppo[min(max(i,2), n-1), j])

# Policy functions
@NLexpression(m, c[i=1:n, j=1:2, s=1:2], abs(Vp[i,j,s])^(-1/γ))

@NLexpression(m, k1[j=1:2, s=1:2], 0.0)
@NLexpression(m, ko[i=2:n, j=1:2, s=1:2], (- Vp[i,j,s] / Vpp[i,j]) * (μR - r) / σR^2)
@NLexpression(m, kk[i=1:n, j=1:2, s=1:2], (i==1) * k1[j,s]  + (i>1) * ko[max(i, 2),j,s])
@NLexpression(m, k[i=1:n, j=1:2, s=1:2], (kk[i,j,s] <= (a[i] - a[1])) * (kk[i,j,s] >= 0) * kk[i, j,s]  + (kk[i,j,s] >= (a[i]-a[1])) * (a[i]-a[1]))

@NLexpression(m, savings[i=1:n, j=1:2, s=1:2],  (z[j] + r * a[i] + (μR - r) * k[i, j, s] - c[i, j, s]))

#Upwind
@NLexpression(m, UVp[i=1:n, j=1:2], (savings[i, j, 1] <= 0) * Vp[i, j, 1] +  (savings[i, j, 1] > 0) * (savings[i, j, 2] >= 0) * Vp[i, j, 2])
@NLexpression(m, Uk[i=1:n, j=1:2], (savings[i, j, 1] <= 0) * k[i, j, 1] +  (savings[i, j, 1] > 0) * (savings[i, j, 2] >= 0) * k[i, j, 2] + (savings[i, j, 1] > 0) * (savings[i, j, 2] < 0) * (k[i,j,1]+k[i,j,2])/2)
@NLexpression(m, Uc[i=1:n, j=1:2], (savings[i, j, 1] <= 0) * c[i, j, 1] +  (savings[i, j, 1] > 0) * (savings[i, j, 2] >= 0) * c[i, j, 2] + (savings[i, j, 1] > 0) * (savings[i, j, 2] < 0) * (z[j] + r * a[i] + (μR-r) * Uk[i,j]))

# Constraint
@NLexpression(m, HJB[i=1:n, j=1:2], Uc[i,j]^(1 - γ) / (1 - γ) + UVp[i,j] * (z[j] + r * a[i] + (μR - r) *Uk[i, j] - Uc[i, j]) + Vpp[i,j] * 0.5 * Uk[i,j]^2 * σR^2 + λ[j] * (V[i, (j==1) * 2 + (j==2) * 1]- V[i, j]))
@NLconstraint(m, CHJB[i=1:n, j=1:2], ρ * V[i, j] == HJB[i, j])



# Objective
@NLobjective(m, Max, 1.0)
@time solve(m)
@time c = getvalue(c)

In REPL mode:
...
julia> @time solve(m)
ERROR: AssertionError: !(nonlinear_wrt_output[nod.parent])
in compute_hessian_sparsity at /home/uma/diam/local/pack/julia/x86_64-darwin/pkgs/v0.4/ReverseDiffSparse/src/sparsity.jl:73
in initialize at /home/uma/diam/local/pack/julia/x86_64-darwin/pkgs/v0.4/JuMP/src/nlp.jl:316
in loadproblem! at /home/uma/diam/local/pack/julia/x86_64-darwin/pkgs/v0.4/Ipopt/src/IpoptSolverInterface.jl:41
in _buildInternalModel_nlp at /home/uma/diam/local/pack/julia/x86_64-darwin/pkgs/v0.4/JuMP/src/nlp.jl:1204
in build at /home/uma/diam/local/pack/julia/x86_64-darwin/pkgs/v0.4/JuMP/src/solvers.jl:307
in solve at /home/uma/diam/local/pack/julia/x86_64-darwin/pkgs/v0.4/JuMP/src/solvers.jl:131

julia> @time c = getvalue(c)

The issue you're having is jump-dev/JuMP.jl#702. Run Pkg.checkout("ReverseDiffSparse")