u[c_, l_] = c^(1 - cgam)/(1 - cgam) + (b*l^(1 - ngam))/(1 - ngam)
 
F[k_, n_, a_] := a*k^alpha*n^(1 - alpha)
 
getSS := Block[{s1, s2, s3, s0}, 
   s1 := {(D[F[k, n, a], k] + 1 - delta/ppy)/(1 + rho/ppy) == 1} /. spec; 
    s2 := {mrs[c, leis] == D[F[k, n, a], n] /. {leis -> 1 - n}} /. spec; 
    s3 := {F[k, n, 1] == c + delta*k/ppy} /. spec; s0 = Join[s1, s2, s3]; 
    sol1 = FindRoot[s0, {c, 0.3}, {k, 1}, {n, 0.3}]; 
    Print["The Steady State has been Found"]; sol1]
 
s1 := {(D[F[k, n, a], k] + 1 - delta/ppy)/(1 + rho/ppy) == 1} /. spec
 
spec = {arho -> 0.95, alpha -> 0.3, cgam -> 1.5, delta -> 0.08, ngam -> 1.5, 
   rho -> 0.02, b -> 2, a -> 1, ppy -> 4}
 
s2 := {mrs[c, leis] == D[F[k, n, a], n] /. {leis -> 1 - n}} /. spec
 
mrs[c_, leis_] := D[u[c, leis], leis]/D[u[c, leis], c] /. spec
 
s3 := {F[k, n, 1] == c + (delta*k)/ppy} /. spec
 
s0 = {0.9950248756218906*(0.98 + (0.3*n^0.7)/k^0.7) == 1, 
   (2*c^1.5)/(1 - n)^1.5 == (0.7*k^0.3)/n^0.3, k^0.3*n^0.7 == c + 0.02*k}
 
sol1 = {c -> 0.6927288482716841, k -> 10.93782384094611, 
   n -> 0.3142269522770712}
 
setFOC := Block[{f1, f2, f3, f4, f5, f6, f7, f8, m0, m1, m2, m3, m4, z1, z2, 
    FK, FN, uc, c2, n2}, Clear[k1, a1]; cbar = c /. sol1; kbar = k /. sol1; 
    nbar = n /. sol1; ibar = k*delta/ppy /. sol1 /. spec; abar = 1; 
    f1 = D[mrs[c, leis], c]*dc + D[mrs[c, leis], leis]*dleis /. 
       {dc -> c1*c, dleis -> -n1*n, leis -> 1 - n} /. {c -> cbar, n -> nbar}; 
    FN[k_, n_, a_] := D[F[k, n, a], n]; 
    f2 = D[FN[k, n, a], k]*dk + D[FN[k, n, a], n]*dn + 
         D[FN[k, n, a], a]*da /. spec /. 
       {dk -> k1*k, dn -> n1*n, da -> a1*a} /. {k -> kbar, n -> nbar}; 
    uc[c_, leis_] := D[u[c, leis], c]; 
    f3 = D[uc[c, leis], c]*dc + D[uc[c, leis], leis]*dleis /. 
        {dc -> c1*c, dleis -> -n1*n, leis -> 1 - n} /. spec /. 
      {c -> cbar, n -> nbar}; FK[k_, n_, a_] := D[F[k, n, a], k]; 
    f4 = Simplify[beta*uc[c, leis]*
            (D[FK[k, n, a], k]*dk + D[FK[k, n, a], n]*dn + 
              D[FK[k, n, a], a]*da) + 
           beta*(D[uc[c, leis], c]*dc + D[uc[c, leis], leis]*dleis)*
            (1 + FK[k, n, a] - delta/ppy) /. {beta -> (1 + rho/ppy)^(-1)} \
/. {dc -> c2*c, dleis -> -n2*n, dk -> k2*k, da -> a2*a, dn -> n2*n, 
          leis -> 1 - n} /. {c -> cbar, n -> nbar, k -> kbar} /. spec]; 
    f5 = k2; f6 = k1*(1 - delta/ppy) + delta/ppy*i1 /. spec; 
    f7 = D[f, k]*dk*k1 + D[f, n]*dn + D[f, a]*da /. 
        {dk -> k1*k, dn -> n1*n, da -> a1*a} /. {k -> kbar, n -> nbar} /. 
      spec; f8 = cbar*c1 + delta/ppy*kbar*i1 /. spec; Clear[m1, m2, m3, m4]; 
    m1 = {f1 == f2}; m2 = {f3 == f4}; m3 = {f5 == f6}; m4 = {f7 == f8}; 
    Clear[c1, n1, i1, k2]; c1 = b11*k1 + b12*a1; n1 := b21*k1 + b22*a1; 
    i1 := b31*k1 + b32*a1; k2 := b41*k1 + b42*a1; c2 = b11*k2 + b12*a2; 
    n2 = b21*k2 + b22*a2; a2 = arho*a1; m0 = Join[m1, m2, m3, m4] /. spec; 
    z1 = Simplify[m0 /. {a1 -> 1, k1 -> 0}]; 
    z2 = Simplify[m0 /. {a1 -> 0, k1 -> 1}]; z0 = Join[z1, z2]; 
    Print["You have set up the FOC's. Ready to solve    for the b's for the \
decision functions."]; ]
 
f1 = 3.045759075433866*c1 + 1.395592309761483*n1
 
c1 = a1*b12 + b11*k1
 
n1 := b21*k1 + b22*a1
 
f2 = 2.03050604814785*a*a1 + 0.6091518144443553*k1 - 0.6091518144443553*n1
 
f3 = -2.60163488971529*c1
 
f4 = 0.04314485741440254*a2 - 2.601634890023081*c2 - 0.03020140019008177*k2 + 
   0.03020140019008177*n2
 
a2 = a1*arho
 
k2 := b41*k1 + b42*a1
 
f5 := k2
 
f6 := k1*(1 - delta/ppy) + delta/ppy*i1 /. spec
 
i1 := b31*k1 + b32*a1
 
f7 := D[F[k, n, a], k]*dk*k1 + D[F[k, n, a], n]*dn + D[F[k, n, a], a]*da /. 
     {dk -> k1*k, dn -> n1*n, da -> a1*a} /. {k -> kbar, n -> nbar} /. spec
 
kbar = 10.93782384094611
 
nbar = 0.3142269522770712
 
f8 := cbar*c1 + delta/ppy*kbar*i1 /. spec
 
cbar = 0.6927288482716841
 
m0 = {3.045759075433866*(a1*b12 + b11*k1) + 
     1.395592309761483*(a1*b22 + b21*k1) == 
    2.03050604814785*a1 + 0.6091518144443553*k1 - 
     0.6091518144443553*(a1*b22 + b21*k1), 
   -2.60163488971529*(a1*b12 + b11*k1) == 
    0.04098761454368241*a1 - 0.03020140019008177*(a1*b42 + b41*k1) - 
     2.601634890023081*(0.95*a1*b12 + b11*(a1*b42 + b41*k1)) + 
     0.03020140019008177*(0.95*a1*b22 + b21*(a1*b42 + b41*k1)), 
   a1*b42 + b41*k1 == 0.98*k1 + 0.02*(a1*b32 + b31*k1), 
   0.9114853244137989*a1 + 0.2734455973241396*k1^2 + 
     0.6380397270896592*(a1*b22 + b21*k1) == 
    0.6927288482716841*(a1*b12 + b11*k1) + 
     0.2187564768189222*(a1*b32 + b31*k1)}
 
m1 := {f1 == f2}
 
m2 := {f3 == f4}
 
m3 := {f5 == f6}
 
m4 := {f7 == f8}
 
z1 = {3.045759075433866*b12 + 1.395592309761483*b22 == 
    2.03050604814785 - 0.6091518144443553*b22, 
   -2.60163488971529*b12 == 0.04098761454368241 - 2.471553145521927*b12 + 
     0.02869133018057769*b22 - 0.03020140019008177*b42 - 
     2.601634890023081*b11*b42 + 0.03020140019008177*b21*b42, 
   b42 == 0.02*b32, 0.9114853244137989 + 0.6380397270896592*b22 == 
    0.6927288482716841*b12 + 0.2187564768189222*b32}
 
z2 = {3.045759075433866*b11 + 1.395592309761483*b21 == 
    0.6091518144443553 - 0.6091518144443553*b21, 
   -2.60163488971529*b11 == (-0.03020140019008177 - 2.601634890023081*b11 + 
       0.03020140019008177*b21)*b41, b41 == 0.98 + 0.02*b31, 
   0.2734455973241396 + 0.6380397270896592*b21 == 
    0.6927288482716841*b11 + 0.2187564768189222*b31}
 
FK[k_, n_, a_] := D[F[k, n, a], k]
 
FN[k_, n_, a_] := D[F[k, n, a], n]
 
uc[c_, leis_] := D[u[c, leis], c]
 
ibar = 0.2187564768189222
 
abar = 1
 
z0 = {3.045759075433866*b12 + 1.395592309761483*b22 == 
    2.03050604814785 - 0.6091518144443553*b22, 
   -2.60163488971529*b12 == 0.04098761454368241 - 2.471553145521927*b12 + 
     0.02869133018057769*b22 - 0.03020140019008177*b42 - 
     2.601634890023081*b11*b42 + 0.03020140019008177*b21*b42, 
   b42 == 0.02*b32, 0 == 0.6927288482716841*b12 + 0.2187564768189222*b32, 
   3.045759075433866*b11 + 1.395592309761483*b21 == 
    0.6091518144443553 - 0.6091518144443553*b21, 
   -2.60163488971529*b11 == (-0.03020140019008177 - 2.601634890023081*b11 + 
       0.03020140019008177*b21)*b41, b41 == 0.98 + 0.02*b31, 
   0 == 0.6927288482716841*b11 + 0.2187564768189222*b31}
 
getbb := Block[{}, sol2 = 
     FindRoot[z0, {b11, 0.2}, {b12, 0.2}, {b21, -0.1}, {b22, 0.2}, 
      {b31, 0.2}, {b32, 0.2}, {b41, 0.2}, {b42, 0.2}]; 
    Print["You have found the b's for the decision      functions. Good."]; 
    sol2]
 
sol2 = {b11 -> 0.326396582557882, b12 -> -0.4913797850567688, 
   b21 -> -0.1920312594628088, b22 -> 1.759391857235889, 
   b31 -> -1.033589185577783, b32 -> 1.55603599727073, 
   b41 -> 0.9593282162884442, b42 -> 0.0311207199454146}
 
preY = Null
 
nextY[x_List] := 
  Block[{k, n, a, y, c, i, newk, newk2, a2, epsilon}, 
   k = kbar*(1 + K1) /. {k1 -> x[[1]], a1 -> x[[2]]}; 
    n = nbar*(1 + N1) /. {k1 -> x[[1]], a1 -> x[[2]]}; 
    c = cbar*(1 + C1) /. {k1 -> x[[1]], a1 -> x[[2]]}; 
    a = abar*(1 + a1) /. {k1 -> x[[1]], a1 -> x[[2]]}; 
    i = ibar*(1 + I1) /. {k1 -> x[[1]], a1 -> x[[2]]}; 
    newk = K2 /. {k1 -> x[[1]], a1 -> x[[2]]}; y = F[k, n, a] /. spec; 
    newk1 = k*(1 - delta/ppy) + i /. spec; newk2 = kbar*(1 + K2); 
    errory = (c + i)/y - 1; errork = (newk1 - newk2)/newk2; 
    epsilon = Random[Real, {-se1, se1}, 2]; 
    a2 = x[[2]]*arho + epsilon /. spec; nextka = {newk, a2}; 
    {k, n, a, y, c, i}]
 
newk2 = 13966.46110296646
 
K1 := k1
 
N1 := b21*k1 + b22*a1 /. sol2
 
C1 = -0.4913797850567688*a1 + 0.326396582557882*k1
 
I1 := b31*k1 + b32*a1 /. sol2
 
K2 := b41*k1 + b42*a1 /. sol2
 
newk1 = 10.93735002148318
 
errory = 0.01078644647367088
 
errork = (0.09142586446277057*(10.93735002148318 - 
       10.93782384094611*(1 + 0.0311207199454146*a1 + 0.9593282162884442*k1))
)/(1 + 0.0311207199454146*a1 + 0.9593282162884442*k1)
 
se1 = 0.002886751345948129
 
nextka = {-0.00004331935399703567, -0.00858630711567012}
 
makeseries[periods_] := 
  Block[{gdp1}, Clear[gdplist]; gdplist = {}; ka = {0, 0}; 
    Do[gdp1 = nextY[ka]; ka = nextka; gdplist = Append[gdplist, gdp1], 
     {periods}]; Print["A matrix called 'gdplist' has been created.  It \
contains in each row the vector     (k,n,a,y,c,i). It is ", periods, 
     " periods long."]; ]
 

 
ka = {-0.00004331935399703567, -0.00858630711567012}