Ok, here it is in all its glory:

```
function f(t, A, prodinput) # have to include parameters as kwargs here
# A is a vector of length 30, prodinput is 160 x 30
HSmale25 = prodinput[1:10, :]
HSmale35 = prodinput[11:20, :]
HSmale45 = prodinput[21:30, :]
HSmale55 = prodinput[31:40, :]
HSfemale25 = prodinput[41:50, :]
HSfemale35 = prodinput[51:60, :]
HSfemale45 = prodinput[61:70, :]
HSfemale55 = prodinput[71:80, :]
COLmale25 = prodinput[81:90, :]
COLmale35 = prodinput[91:100, :]
COLmale45 = prodinput[101:110, :]
COLmale55 = prodinput[111:120, :]
COLfemale25 = prodinput[121:130, :]
COLfemale35 = prodinput[131:140, :]
COLfemale45 = prodinput[141:150, :]
COLfemale55 = prodinput[151:160, :]
for i in 1:3, k in 1:2 # skilled: occupations 1,2
τ[i,k,1] = τ[i,1,1]
σ[i,k,1] = σ[i,1,1]
ν[i,k,1] = ν[i,1,1]
ϖ[i,k,1] = ϖ[i,1,1]
end
for i in 1:3, k in 3:6 # services: occupations 3-6
τ[i,k,1] = τ[i,1,1]
σ[i,k,1] = σ[i,1,1]
ν[i,k,1] = ν[i,1,1]
ϖ[i,k,1] = ϖ[i,1,1]
end
for i in 1:3, k in 7:10 # services: occupations 7-10
τ[i,k,1] = τ[i,1,1]
σ[i,k,1] = σ[i,1,1]
ν[i,k,1] = ν[i,1,1]
ϖ[i,k,1] = ϖ[i,1,1]
end
for i in 1:3, k in 1:10, j in 1:5
τ[i,k,j] = τ[i,1,j]
σ[i,k,j] = σ[i,1,j]
ν[i,k,j] = ν[i,1,j]
ϖ[i,k,j] = ϖ[i,1,j]
end
for t in 1:szT
# equations 8a-8d
for k in 1:szK
HSmale[k, t] = (logit(τ[1,k,:],t)*HSmale25[k, t]^ρ[7] + logit(τ[2,k,:],t)*HSmale35[k, t]^ρ[7] + logit(τ[3,k,:],t)*HSmale45[k, t]^ρ[7] + (1-logit(τ[1,k,:],t)-logit(τ[2,k,:],t)-logit(τ[3,k,:],t))*HSmale55[k, t])^(1/ρ[7])
HSfemale[k, t] = (logit(σ[1,k,:],t)*HSfemale25[k, t]^ρ[7] + logit(σ[2,k,:],t)*HSfemale35[k, t]^ρ[7] + logit(σ[3,k,:],t)*HSfemale45[k, t]^ρ[7] + (1-logit(σ[1,k,:],t)-logit(σ[2,k,:],t)-logit(σ[3,k,:],t))*HSfemale55[k, t])^(1/ρ[7])
COLmale[k, t] = (logit(ν[1,k,:],t)*COLmale25[k, t]^ρ[7] + logit(ν[2,k,:],t)*COLmale35[k, t]^ρ[7] + logit(ν[3,k,:],t)*COLmale45[k, t]^ρ[7] + (1-logit(ν[1,k,:],t)-logit(ν[2,k,:],t)-logit(ν[3,k,:],t))*COLmale55[k, t])^(1/ρ[7])
COLfemale[k, t] = (logit(ϖ[1,k,:],t)*COLfemale25[k, t]^ρ[7] + logit(ϖ[2,k,:],t)*COLfemale35[k, t]^ρ[7] + logit(ϖ[3,k,:],t)*COLfemale45[k, t]^ρ[7] + (1-logit(ϖ[1,k,:],t)-logit(ϖ[2,k,:],t)-logit(ϖ[3,k,:],t))*COLfemale55[k, t])^(1/ρ[7])
end
for k in 1:szK # eqs 7
HS[k, t] = (logit(ξ[k,:],t)*HSmale[k, t]^ρ[6] + (1-logit(ξ[k,:],t))*HSfemale[k, t]^ρ[6])^(1/ρ[6])
COL[k, t] = (logit(γ[k,:],t)*COLmale[k, t]^ρ[6] + (1-logit(γ[k,:],t))*COLfemale[k, t]^ρ[6])^(1/ρ[6])
end
for k in 1:szK # eq 6
EMPk[k, t] = (logit(μ[k,:],t)*HS[k, t]^ρ[5] + (1-logit(μ[k,:],t))*COL[k, t]^ρ[5])^(1/ρ[5])
end
# eqs 5
EMPu[1, t] = (logit(λ[3,:],t)*EMPk[3, t]^ρ[3] + logit(λ[4,:],t)*EMPk[4, t]^ρ[3] + logit(λ[5,:],t)*EMPk[5, t]^ρ[3] + (1-logit(λ[3,:],t)-logit(λ[4,:],t)-logit(λ[5,:],t))*EMPk[6, t]^ρ[3])^(1/ρ[3])
EMPu[2, t] = (logit(λ[7,:],t)*EMPk[7, t]^ρ[4] + logit(λ[8,:],t)*EMPk[8, t]^ρ[4] + logit(λ[9,:],t)*EMPk[9, t]^ρ[4] + (1-logit(λ[7,:],t)-logit(λ[8,:],t)-logit(λ[9,:],t))*EMPk[10, t]^ρ[4])^(1/ρ[4])
# eqs 4
EMP_su[1, t] = (logit(λ[1,:],t)*EMPk[1, t]^ρ[1] + (1-logit(λ[1,:],t))*EMPk[2, t]^ρ[1])^(1/ρ[1])
EMP_su[2, t] = (logit(λ[2,:],t)*EMPu[1, t]^ρ[2] + (1-logit(λ[2,:],t))*EMPu[2, t]^ρ[2])^(1/ρ[2])
end
Y = β(t; b=b)*(λY[2]*(λY[1]*A[t]^ρY[1] + (1 - λY[1])*EMP_su[1,t]^ρY[1])^(ρY[2]/ρY[1]) + (1-λY[2])*EMP_su[2,t]^ρY[2])^(1/ρY[2])
return Y
end
```

Now that I look at it after giving it a facelift, I think it can work with `ForwardDiff`

to obtain the derivatives wrt `prodinput`

, no?