General Setup

Needs["Statistics`NormalDistribution`"]
Off[General::spell]
Off[General::spell1]

Calibration

Clear @@ 
      Hold[α,  β,  δ,  θ, τ, ρ, σ]

α = 0.4;
β = 0.9896;
δ = 0.0196;
θ = 0.357;
τ = 2;
ρ = 0.95;
σ = 0.007;

Steady State

ϕ = (((((((1/β) - 1 + δ))/α)))^((1/(1 - α))));
ω = (ϕ^(1 - α)) - δ;
Φ = (θ/(1 - θ)) ((1 - α)/(ϕ^α));
kss = (Φ/(ω + ϕ Φ));
lss = ϕ kss;
css = ω kss;
yss = (kss^α) 
          (lss^(1 - α));
(Subscript[α, 1]) = ((θ (1 - τ) - 1)/css);
(Subscript[α, 2]) = -(((1 - θ) (1 - τ))/(1 - 
                lss));
(Subscript[α, 3]) = β (( α (1 - α))/lss) (kss^(α - 1)) 
            (lss^(1 - α)) + (Subscript[α, 2]);
(Subscript[α, 4]) = α β (kss^(α - 1)) 
          (lss^(1 - α));
(Subscript[α, 5]) = β (( α (α - 
                    1))/kss) (kss^(α - 1)) 
          (lss^(1 - α));
(Subscript[α, 6]) = -((α/lss) + (1/(1 - lss))) 
          css;
(Subscript[α, 7]) = 
        yss (( 1 - α)/lss);
A = 1;
B = (α/kss) (css - 
                yss) - (1 - δ);
CC = (Subscript[α, 6]) - (Subscript[α, 7]);
DD = css - yss;
G = (Subscript[α, 1]) (α/kss) 
            css + (Subscript[α, 5]);
H = -(Subscript[α, 1]) (α/kss) 
          css;
J = (Subscript[α, 1]) (Subscript[α, 6]) + (Subscript[α, 3]);
K = -((Subscript[α, 1]) (Subscript[α, 6]) + (Subscript[α, 2]));
L = (Subscript[α, 1]) css + (Subscript[α, 4]);
M = -(Subscript[α, 1]) css;
NN = ρ;
sol = Solve[{(x^2) + ((B/A) + (K/J) - ((G CC)/(J A))) 
                    x + ((K B  - H CC)/(J A)) == 
                0}, {x}][[1]];
P = x  /.  sol;
R = -(1/CC) (A P + B);
Q = ((-DD (J NN  + K) + CC L NN  + CC M)/(A J NN + A K - CC G - CC J R));
S = ((-A L NN - A M + DD G  + DD J R)/(A J NN + A K - CC G - CC J R));

sslist = {kss, lss, P, Q, R, S}

Tauchen and Capital Grid

Tauchen

cdf[x_] := CDF[NormalDistribution[0, 1], x]
points = 5;

    
\(\(\(Tpoints[a_, points_] := 
      Table[\(-3\)*a + 2*3*a*\((n - 1)\)/\((points - 1)\), {n, 
          points}]\)\(
\)
    \)\)
Tmatrix[A_, σ_, ρ_, 
        points_] := 
Table[
Which[
j == 
            1, 
cdf[( 
                1/2*( 
                      A[[j + 1]] + 
                        A[[j]]) - ρ A[[i]])/σ], 
j == points, 
1 - 
            cdf[( 1/2*( 
                        A[[j]] + 
                          A[[j - 
                                1]]) - ρ A[[i]])/σ], 
True, 
cdf[( 
                  1/2*( 
                        A[[j + 1]] + 
                          A[[j]]) - ρ A[[i]])/σ] - 
            cdf[( 1/2*( 
                        A[[j - 1]] + 
                          A[[j]]) - ρ A[[i]])/σ]
], 
{i, points}, {j, 
          points}]

Capital Grid

coverk = 0.299;
capitalpoints = 5;
CP = 
        Table[kss*(1 - coverk) + (i - 1)*kss*2*
              coverk/(capitalpoints - 1), {i, 
            capitalpoints}];

Functions

KapitalTomorrow[k_, 
          z_, γ_, ζ_, φ_, {kss_, lss_, P_, Q_, R_, 
            S_}] := (((((γ/ζ) (kss^(γ - ζ)) P ((k^ζ) - (kss^ζ)) + (γ/φ) (kss^(γ - 1)) Q*(z^φ)) + (kss^γ)))^((1/γ)))

Labor[k_, 
          z_, μ_, ζ_, φ_, {kss_, lss_, P_, Q_, R_, 
            S_}] := (((((μ/ζ) (lss^(μ - 1)) (kss^(1 - ζ)) R ((k^ζ) - (kss^ζ)) + (μ/φ) (lss^(μ - 1)) S*(z^φ)) + (lss^μ)))^((1/μ)))

Output[k_, 
          z_, μ_, ζ_, φ_, {kss_, lss_, P_, Q_, R_, 
            S_}] := (kss^α) 
            (lss^(1 - α)) + (kss^α) 
            (lss^(1 - α))*(z + α*(k - kss)/
                    kss + (1 - α)*(Labor[k, 
                          z, μ, ζ, φ, {kss, lss, P, Q, R, 
                            S}] - lss)/lss)

Consumption[k_,z_,γ_,μ_,ζ_,φ_,{kss_, lss_, P_, Q_, R_, S_}]:=Output[k, 
            z, μ, ζ, φ, {kss, lss, P, Q, R, 
              S}] + (1 - δ)*k - 
          KapitalTomorrow[k, 
            z, γ, ζ, φ, {kss, lss, P, Q, R, S}]

InteresRate[k_, 
          z_, μ_, ζ_, φ_, {kss_, lss_, P_, Q_, R_, S_}] := 
        1 - δ + α*Exp[z]*(k^(α - 1)) 
            Labor[k, z, μ, ζ, φ, {kss, lss, P, Q, R, S}]^(1 - α)

InterestRatetomorrow[k_, z_, 
          zp_, γ_, μ_, ζ_, φ_, {kss_, lss_, P_, 
            Q_, R_, S_}] := 
        InteresRate[
          KapitalTomorrow[k, 
            z, γ, ζ, φ, {kss, lss, P, Q, R, S}], 
          zp, μ, ζ, φ, {kss, lss, P, Q, R, S}]

Utility[c_, 
          l_] := (((((c^θ) (((1 - l))^(1 - θ))))^(1 - τ))/(1 - τ))

MarginalUtilityConsumption[c_, 
          l_] := ((θ (1 - τ) Utility[c, l])/c)

TomorrowMarginalUtilityInterestRate[k_, z_, 
          zp_, γ_, μ_, ζ_, φ_, {kss_, lss_, P_, 
            Q_, R_, S_}]:=β*MarginalUtilityConsumption[Consumption[KapitalTomorrow[k, 
                    z, γ, ζ, φ, {kss, lss, P, Q, R, 
                      S}],zp,γ,μ,ζ,φ,{kss, lss, P, Q, R, S}],Labor[
                KapitalTomorrow[k, 
                  z, γ, ζ, φ, {kss, lss, P, Q, R, S}], 
                zp, μ, ζ, φ, {kss, lss, P, Q, R, S}]]*InteresRate[
            KapitalTomorrow[k, 
              z, γ, ζ, φ, {kss, lss, P, Q, R, S}], 
            zp, μ, ζ, φ, {kss, lss, P, Q, R, S}]

TomorrowMarginalUtilityInterestRateExpectation[k_, TP_, TM_, 
        i_, γ_, μ_, ζ_, φ_, {kss_, lss_, P_, Q_, 
          R_, S_}] := (∑_(j = 1)^(Length[TP]))TM[[i, j]] 
          TomorrowMarginalUtilityInterestRate[k, TP[[i]], 
            TP[[j]], γ, μ, ζ, φ, {kss, lss, P, 
              Q, R, S}]
InverseUtility[A_, c_, 
          l_] := (((((A/(θ (((1 - l))^((1 - θ) (1 - τ)))))))^((1/(θ (1 - τ) - 1))))/c)

EulerErrors[k_, TP_, TM_, 
        i_, γ_, μ_, ζ_, φ_, {kss_, lss_, P_, Q_, 
          R_, S_}]:=1-InverseUtility[TomorrowMarginalUtilityInterestRateExpectation[k, TP, TM, 
              i, γ, μ, ζ, φ, {kss, lss, P, Q, R, 
                S}],Consumption[k,TP[[i]],γ,μ,ζ,φ,{kss, lss, P, Q, R, S}],Labor[k, 
              TP[[i]], μ, ζ, φ, {kss, lss, P, Q, R, 
                S}]]

Results

Euler[γ_, μ_, ζ_, φ_, {kss_, lss_, P_, Q_, 
          R_, S_}] := 
      Table[Abs[
          EulerErrors[CP[[i]], TP, TM, 
            j, γ, μ, ζ, φ, {kss, lss, P, Q, R, 
              S}]], {i, 1, capitalpoints}, {j, 1, points}]

TotalEuler[γ_, μ_, ζ_, φ_, {kss_, 
            lss_, P_, Q_, R_, S_}] := 
        Total[Euler[γ, μ, ζ, φ, {kss, lss, P, Q, 
              R, S}], 3]

Minimization with 21 (points and capitalpoints)

points = 21;
capitalpoints = 21;
TP = Tpoints[σ/(Sqrt[1 - (ρ^2)]), points];
TM = 
        Tmatrix[TP, σ, ρ, points] // 
          Chop;
coverk = 0.299;
CP = 
        Table[kss*(1 - coverk) + (i - 1)*kss*2*
              coverk/(capitalpoints - 1), {i, 
            capitalpoints}];

TotalEuler[1, 1, 1, 1, sslist] // Timing
list = {};
fm1  =  
      FindMinimum[
          TotalEuler[γ, μ, ζ, 1, sslist], {γ, 
            0.9587004360294477}, {ζ, 0.9673663282097958}, {μ, 
            2.477985230984596}, 
          StepMonitor :> (Print[
                TotalEuler[γ, μ, ζ, 1, sslist]]; 
              AppendTo[list, {γ, μ}]), Compiled -> False] // 
        Timing