diff --git a/Project.toml b/Project.toml
index 26a66f6..7d78260 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Florian Oswald <florian.oswald@gmail.com>"]
 version = "0.1.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
diff --git a/runb.fish b/runb.fish
new file mode 100755
index 0000000..a5e6d69
--- /dev/null
+++ b/runb.fish
@@ -0,0 +1,4 @@
+ #! /usr/local/bin/fish
+for n in (seq 4)
+    env JULIA_NUM_THREADS=$n julia --project=. ./src/bmark.jl   
+end
diff --git a/src/bmark.jl b/src/bmark.jl
new file mode 100644
index 0000000..2b7374d
--- /dev/null
+++ b/src/bmark.jl
@@ -0,0 +1,14 @@
+using DCEGM
+using BenchmarkTools
+println("using $(Threads.nthreads()) threads")
+println("===========================")
+na = 5000; ny = 50
+
+println("fedors model on na=$na,ny=$ny")
+DCEGM.runf(par = Dict(:na => na, :ny => ny));
+@btime DCEGM.runf(par = Dict(:na => na, :ny => ny));
+
+na = 500; ny = 15
+println("general model on na=$na,ny=$ny")
+DCEGM.rung(par = Dict(:na => na, :ny => ny));
+@btime DCEGM.rung(par = Dict(:na => na, :ny => ny));
diff --git a/src/dc_algo.jl b/src/dc_algo.jl
index ad0be96..f23469e 100644
--- a/src/dc_algo.jl
+++ b/src/dc_algo.jl
@@ -110,6 +110,7 @@ function dc_EGM!(m::FModel,p::Param)
             end
         else
             for id in 1:p.nD   # current period dchoice
+            # Threads.@threads for id in 1:p.nD   # current period dchoice
                 working = id==1  # working today is id=1
                 # what's next period's cash on hand given you work/not TODAY?
                 mm1 = m.m1[it][id]
@@ -121,7 +122,7 @@ function dc_EGM!(m::FModel,p::Param)
                 vmat = fill(-Inf,p.nD,p.na*p.ny)
 
                 # get future cons and values
-                for iid in 1:p.nD
+                Threads.@threads for iid in 1:p.nD
                     c1 = interp(m.c[iid,it+1].env, mm1[:])
                     floory!(c1,p.cfloor)   # floor negative consumption
                     cmat[iid,:] = gety(c1)  # get y-values
@@ -191,7 +192,7 @@ function dc_EGM!(m::GModel,p::Param)
     cmat = fill(-Inf,p.nD,p.na)
     vmat = fill(-Inf,p.nD,p.na)
     ctmp = fill(-Inf,p.nD,p.ny,p.na)
-    # vtmp = fill(-Inf,p.nD,p.ny,p.na)
+    vtmp = fill(-Inf,p.ny,p.nD,p.na)
 
     for it in p.nT:-1:1
         for iy in 1:p.ny  # current state
@@ -210,12 +211,14 @@ function dc_EGM!(m::GModel,p::Param)
                     # reset all value matrices
                     fill!(vmat,0.0)
                     fill!(ctmp,-Inf)
-                    # fill!(vtmp,0.0)
+                    # fill!(vtmp,-Inf)
 
+                    # Threads.@threads for jy in 1:p.ny # future state: owner, renter, income, etc
                     for jy in 1:p.ny # future state: owner, renter, income, etc
                         pr = m.ywgt[iy,jy]  # transprob
 
-                        for iid in 1:p.nD  # future dchoice
+                        Threads.@threads for iid in 1:p.nD  # future dchoice
+                        # for iid in 1:p.nD  # future dchoice
                             # only feasible choices at this state
                             # if renter, cannot sell etc
 
@@ -223,14 +226,17 @@ function dc_EGM!(m::GModel,p::Param)
                             c1 = interp(m.c[iid,jy,it+1].env, m1) # C(d',y',m')
                             floory!(c1,p.cfloor)   # floor negative consumption
                             ctmp[iid,jy,:] = gety(c1)
-                            # vtmp[iid,jy,:] = vfun(iid,it+1,ctmp[iid,jy,:],m1,m.v[iid,jy,it+1],p)
+                            # vtmp[jy,iid,:] = pr * vfun(iid,it+1,ctmp[iid,jy,:],m1,m.v[iid,jy,it+1],p)
                             vmat[iid,:] += pr * vfun(iid,it+1,ctmp[iid,jy,:],m1,m.v[iid,jy,it+1],p)
                         end
                     end # end future state
 
                     # now get expectated value function conditional on iy: E[V(t+1,iid,y')|iy]
-                    # vmat = integrate over second dimension
-                    # vmat = dropdims( reduce(+, vtmp, dims = 2), dims = 2)
+                    # vmat = integrate over first dimension
+                    # vmat2 = dropdims( reduce(+, vtmp, dims = 1), dims = 1)
+                    # @assert all(vmat2 .== vmat)
+
+                    # vmat = dropdims( reduce(+, vtmp, dims = 1), dims = 1)
 
                     # get ccp of choices: P(d'|iy), pwork
                     pwork = working ? ccp(vmat,p) : zeros(size(vmat)[2])