asymptotics work together with winding criterion and give the exact r…

…ight answer, also "saved" is in workflow now

determinacy.py

@@ -0,0 +1,60 @@
+import numpy as np
+from numpy.fft import fftn, rfftn
+import jacobian
+from numba import njit
+
+
+def detA_path(A, N=4096):
+    """Same as detA_path_simple but uses conjugate symmetry across real axis to halve work"""
+    # preliminary: assume and verify shape 2*T-1, k, k for A
+    T = (A.shape[0]+1) // 2
+    k = A.shape[1]
+    assert T == (A.shape[0]+1)/2 and N >= 2*T-1 and k == A.shape[2]
+
+    # step 1: use FFT to calculate A(lambda) for each lambda = 2*pi*{0, 1/N, ..., 1/2} (last if N even)
+    Alambda = rfftn(A[::-1,...], axes=(0,), s=(N,))
+
+    # step 2: take determinant of each and adjust for T-1 displacement
+    det_Alambda = np.empty(N+1, dtype=np.complex128)
+    det_Alambda[:N//2+1] = np.linalg.det(Alambda)*np.exp(2j*np.pi*k*(T-1)/N*np.arange(N//2+1))
+
+    # step 3: use conjugate symmetry to fill in rest (only not duplicating 1/2 if N even)
+    det_Alambda[N//2+1:] = det_Alambda[:(N+1)//2][::-1].conj()
+
+    return det_Alambda
+
+
+@njit
+def winding_number(x, y):
+    """compute winding number of (x,y) coordinates that make closed path by counting
+    number of counterclockwise crossings of ray from (0,0) -> (infty,0) on x axis"""
+    # ensure closed path!
+    assert x[-1] == x[0] and y[-1] == y[0]
+
+    winding_number = 0
+    cur_sign = (y[0] >= 0)
+    for i in range(1, len(x)):
+        if (y[i] >= 0) != cur_sign:
+            # we only enter this if statement rarely (when x axis crossed)
+            # so efficiency no biggie
+            cur_sign = (y[i] >= 0)
+
+            # possible crossing, let's test the three cases
+            if x[i] > 0 and x[i-1] > 0:
+                # entirely on right half-plane, definite crossing
+                winding_number += 2*cur_sign-1
+            elif not (x[i] <= 0 and x[i-1] <= 0):
+                # if not entirely on left half-plane, ambiguous, must check criterion
+                # this step is intended to be rare
+                cross_coord = (x[i-1]*y[i] - x[i]*y[i-1])/(y[i]-y[i-1])
+                if cross_coord > 0:
+                    winding_number += 2*cur_sign-1
+    return winding_number
+
+
+def winding_criterion(A, N=4096):
+    """Build path of det A(lambda) and obtain its winding number, implementing extended
+    Onatski criterion for determinacy: 0 for determinate solution, -1 (or lower)
+    for indeterminacy, 1 (or higher) for no solution"""
+    det_Alambda = detA_path(A, N)
+    return winding_number(det_Alambda.real, det_Alambda.imag)

het_block.py

@@ -326,8 +326,10 @@ def ajac(self, ss, T, shock_list, output_list=None, h=1E-4, Tpost=None, save=Fal
         if use_saved and 'curlyYs' not in self.saved:
             asympJ = {}
             for o in output_list:
+                asympJ[o.upper()] = {}
                 for i in shock_list:
-                    asympJ[o.upper()][i] = self.saved['asympJ'][o.upper()][i][-(Tpost-1): Tpost]
+                    asympJ[o.upper()][i] = asymptotic.AsymptoticTimeInvariant(
+                                                self.saved['asympJ'][o.upper()][i][-(Tpost-1): Tpost])
             return asympJ
 
         # was either saved last by jac or not saved at all, need to do more work!

jacobian.py

@@ -16,7 +16,7 @@
 '''
 
 
-def get_H_U(block_list, unknowns, targets, T, ss=None, asymptotic=False, Tpost=None):
+def get_H_U(block_list, unknowns, targets, T, ss=None, asymptotic=False, Tpost=None, save=False, use_saved=False):
     """Get n_u*T by n_u*T matrix H_U, Jacobian mapping all unknowns to all targets.
 
     Parameters
@@ -33,7 +33,7 @@ def get_H_U(block_list, unknowns, targets, T, ss=None, asymptotic=False, Tpost=N
     """
 
     # do topological sort and get curlyJs
-    curlyJs, required = curlyJ_sorted(block_list, unknowns, ss, T, asymptotic, Tpost)
+    curlyJs, required = curlyJ_sorted(block_list, unknowns, ss, T, asymptotic, Tpost, save, use_saved)
 
     # do matrix forward accumulation to get H_U = J^(curlyH, curlyU)
     H_U_unpacked = forward_accumulate(curlyJs, unknowns, targets, required)
@@ -48,7 +48,7 @@ def get_H_U(block_list, unknowns, targets, T, ss=None, asymptotic=False, Tpost=N
 
 
 def get_impulse(block_list, dZ, unknowns, targets, T=None, ss=None, outputs=None,
-                                                            H_U=None, H_U_factored=None):
+                                        H_U=None, H_U_factored=None, save=False, use_saved=False):
     """Get a single general equilibrium impulse response.
 
     Extremely fast when H_U_factored = utils.factor(get_HU(...)) has already been computed
@@ -76,7 +76,8 @@ def get_impulse(block_list, dZ, unknowns, targets, T=None, ss=None, outputs=None
             T = len(x)
             break
 
-    curlyJs, required = curlyJ_sorted(block_list, unknowns + list(dZ.keys()), ss, T)
+    curlyJs, required = curlyJ_sorted(block_list, unknowns + list(dZ.keys()), ss, T,
+                                      save=save, use_saved=use_saved)
 
     # step 1: if not provided, do (matrix) forward accumulation to get H_U = J^(curlyH, curlyU)
     if H_U is None and H_U_factored is None:
@@ -112,7 +113,7 @@ def get_impulse(block_list, dZ, unknowns, targets, T=None, ss=None, outputs=None
 
 
 def get_G(block_list, exogenous, unknowns, targets, T, ss=None, outputs=None, 
-                                                        H_U=None, H_U_factored=None):
+                            H_U=None, H_U_factored=None, save=False, use_saved=False):
     """Compute Jacobians G that fully characterize general equilibrium outputs in response
     to all exogenous shocks in 'exogenous'
 
@@ -139,7 +140,8 @@ def get_G(block_list, exogenous, unknowns, targets, T, ss=None, outputs=None,
     """
 
     # step 1: do topological sort and get curlyJs
-    curlyJs, required = curlyJ_sorted(block_list, unknowns + exogenous, ss, T)
+    curlyJs, required = curlyJ_sorted(block_list, unknowns + exogenous, ss, T,
+                                      save=save, use_saved=use_saved)
 
     # step 2: do (matrix) forward accumulation to get
     # H_U = J^(curlyH, curlyU) [if not provided], H_Z = J^(curlyH, curlyZ)
@@ -165,7 +167,7 @@ def get_G(block_list, exogenous, unknowns, targets, T, ss=None, outputs=None,
     return forward_accumulate(curlyJs, exogenous, outputs, required | set(unknowns))
 
 
-def curlyJ_sorted(block_list, inputs, ss=None, T=None, asymptotic=False, Tpost=None):
+def curlyJ_sorted(block_list, inputs, ss=None, T=None, asymptotic=False, Tpost=None, save=False, use_saved=False):
     """
     Sort blocks along DAG and calculate their Jacobians (if not already provided) with respect to inputs
     and with respect to outputs of other blocks
@@ -195,9 +197,11 @@ def curlyJ_sorted(block_list, inputs, ss=None, T=None, asymptotic=False, Tpost=N
             jac = block.jac(ss, shock_list=[i for i in block.inputs if i in shocks])
         elif isinstance(block, het.HetBlock):
             if asymptotic:
-                jac = block.ajac(ss, T=T, shock_list=[i for i in block.inputs if i in shocks], Tpost=Tpost)
+                jac = block.ajac(ss, T=T,
+                                 shock_list=[i for i in block.inputs if i in shocks], Tpost=Tpost, save=save, use_saved=use_saved)
             else:
-                jac = block.jac(ss, T=T, shock_list=[i for i in block.inputs if i in shocks])
+                jac = block.jac(ss, T=T,
+                                shock_list=[i for i in block.inputs if i in shocks], save=save, use_saved=use_saved)
         else:
             jac = block
         curlyJs.append(jac)

simple_block.py

@@ -260,7 +260,12 @@ def __add__(self, A):
             return A.reshape((T, T))
 
     def __radd__(self, A):
-        return self + A
+        try:
+            return self + A
+        except:
+            print(self)
+            print(A)
+            raise
 
     def __sub__(self, A):
         # slightly inefficient implementation with temporary for simplicity