WIP

Author: albop

ICAI/partial.dyno ICAI/partial_deterministic.dynoICAI/partial.dyno renamed to ICAI/partial_deterministic.dyno

@@ -51,13 +51,14 @@ c_τ[~] <- b_τ[~]*r + y_τ[~]
 
 
 # exogenous processes
+
 y[t] = y[~] + ρ_y*(y[t-1]-y[~])
 ζ[t] = ζ[~] + ρ_ζ*(ζ[t-1]-ζ[~])
 
 a[t] = 1 + ρ_a*(a[t-1]-1)
 
-d[t]   = d[t-1]   + e_d[t]    # income shock for top earners (aggregate)
-d_b[t] = d_b[t-1] + e_d_b[t]  # income shock for bottom earners (aggregate)
+d[t]   = 0.99999*d[t-1]   + e_d[t]    # income shock for top earners (aggregate)
+d_b[t] = 0.99999*d_b[t-1] + e_d_b[t]  # income shock for bottom earners (aggregate)
 
 y_b[t] = a[t]*(1-y[t])*(1-ζ[t])                      # labor income for bottom-earners (aggregate)
 y_τ[t] = a[t]*(1-y[t])*ζ[t] + (a[t]*y[t])            # top earners income (per capita)
@@ -66,7 +67,7 @@ c_b[t] = d_b[t] +  b_b[t-1]*r - (b_b[t]-b_b[t-1])*p[t] + y_b[t]  # consumption b
 c_τ[t] =   d[t] +  b_τ[t-1]*r - (b_τ[t]-b_τ[t-1])*p[t] + y_τ[t]  # consumption top earners    (aggregate)
 
 W_τ[t] = q[t] + p[t]*(b_τ[t])     # wealth top earners (aggregate)
-W_b[t] = 1         + p[t]*(b_b[t])     # wealth bottom earners (aggregate)
+W_b[t] = 1    + p[t]*(b_b[t])     # wealth bottom earners (aggregate)
 
 p[t] = β_b*(c_b[t+1]/c_b[t])^(-σ)*(r+p[t+1]) + φ_b*(W_b[t])^(-η_b)/(c_b[t])^(-σ)*p[t]
 p[t] = β_τ*(c_τ[t+1]/c_τ[t])^(-σ)*(r+p[t+1]) + φ_τ*(W_τ[t])^(-η  )/(c_τ[t])^(-σ)*p[t]
@@ -83,11 +84,12 @@ r_w + p[t] = 1
 
 
 ### Exogenous shocks
-# e_d[t]   <- N(0, 0.01)
-# e_d_b[t] <- N(0, 0.01)
+# e_d[t]   <- N(0, 0.01^2)
+# e_d_b[t] <- N(0, 0.01^2)
 
 # exogenous shocks
 e_d[1] <- 0.01
+e_d[2] <- 0.0
 e_d_b[1] <- 0.0 # dummy assignment
 
 # forall t, 1<=t<T+1: e_d[t] <- 0.01

ICAI/solve_icai.py

@@ -4,7 +4,7 @@
 from dyno.symbolic_model import SymbolicModel
 import pandas as pd
 
-model = SymbolicModel("ICAI/partial.dyno")
+model = SymbolicModel("ICAI/partial_deterministic.dyno")
 model2  = SymbolicModel("ICAI/partial_shock.dyno")
 dr2 = model2.solve()
 sim2 = dr2.irfs(type='level')['e_d'] # level seems to be ignored
@@ -17,13 +17,28 @@
 pd.DataFrame(v0, columns=model.symbols['variables'])
 
 from dyno.solver import deterministic_solve
+from dyno.solver import solve as linsolve
 import time
 
 
-deterministic_solve(model)
+r0,A,B,C,D = model.jacobians
+X, evs = linsolve(A,B,C)
+r, J = model.deterministic_residuals_with_jacobian(v0.ravel(), sparsify=True)
+rr = r.reshape(v0.shape)
+
+
+
+
+
+import numpy as np
+
+np.isnan(r).sum()
+np.isnan(J.todense()).sum()
+
+from scipy.sparse.linalg import spsolve
 
 sims = [sim2]
-for t in [100,150,200]:
+for t in [150,200,250]:
 
     m = model.recalibrate(T=t)
 
@@ -33,27 +48,44 @@
     print("Elapsed time: ", t2 - t1)
     sims.append(sim)
 
+sims.reverse()
+
+
+
+
 
 plt.figure()
 
 for i,sim in enumerate(sims):
-    plt.subplot(131)
+    plt.subplot(221)
     plt.plot(sim.index, sim['W_τ'], label=i)
+    plt.ylim(2.4, 2.5)
     plt.legend()
     plt.xlim(0,50)
-    plt.subplot(132)
+    plt.subplot(222)
     plt.plot(sim.index, sim['c_τ'])
     plt.xlim(0,50)
-    plt.subplot(133)
+    plt.subplot(223)
     plt.plot(sim.index, sim['d'])
-    plt.xlim(0,50)
+    try:
+        plt.subplot(224)
+        plt.plot(sim.index, sim['e_d'])
+        plt.ylim(-0.001, 0.01)
+    except:
+        pass
+
 
+plt.xlim(0,20)
 plt.tight_layout()
 
 plt.legend()
 
 
 
+
+
+
+
 model.jacobians[0]
 
 

dev_import_tests.py

@@ -4,11 +4,8 @@
 
 from dyno.symbolic_model import SymbolicModel
 
-
 data = yaml.safe_load(open("import_tests.yaml"))
-
 model = SymbolicModel(txt=data['solvable']['rbc_deterministic'])
-
 model
 
 

src/dyno/model.py

@@ -227,6 +227,7 @@ def solve(self: Self, **args):
             return sol
         else:
             dr = self.perturb(**args)
+            return dr
 
     def perturb(self: Self, method: Solver = "qz") -> RecursiveSolution:
         """linearizes the model

src/dyno/solver.py

@@ -499,7 +499,8 @@ def deterministic_solve(model, x0=None, T=None, method="hybr", verbose=True, **a
         u0,
         jactype="sparse",
         verbose=verbose,
-        maxit = args.get("maxit",10)
+        maxit = args.get("maxit",10),
+        tol=args.get("tol",1e-8)
     )
 
     w0 = res.reshape(v0.shape)