Skip to content

Commit

Permalink
commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@marcdelabarrera
marcdelabarrera committed Dec 1, 2022
1 parent a12c783 commit 070a8aa
Show file tree
Hide file tree
Showing 9 changed files with 523 additions and 45 deletions.
45 changes: 40 additions & 5 deletions notebooks/dsolve_examples.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -34,7 +34,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 2,
"metadata": {},
"outputs": [
{
Expand All @@ -56,7 +56,7 @@
"E_{t}[\\pi^{w}_{t+1}]"
]
},
"execution_count": 10,
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -82,7 +82,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 3,
"metadata": {},
"outputs": [
{
Expand All @@ -94,7 +94,7 @@
" 0.13421773, 0.10737418])}"
]
},
"execution_count": 11,
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
Expand Down Expand Up @@ -216,6 +216,41 @@
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Constants"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"ename": "KeyError",
"evalue": "'A_0'",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mKeyError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32mc:\\Users\\MBBar\\OneDrive\\Escritorio\\git_local\\dsolve\\notebooks\\dsolve_examples.ipynb Cell 12\u001b[0m in \u001b[0;36m<cell line: 5>\u001b[1;34m()\u001b[0m\n\u001b[0;32m <a href='vscode-notebook-cell:/c%3A/Users/MBBar/OneDrive/Escritorio/git_local/dsolve/notebooks/dsolve_examples.ipynb#X20sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m eq \u001b[39m=\u001b[39m [\u001b[39m'\u001b[39m\u001b[39mx_\u001b[39m\u001b[39m{t}\u001b[39;00m\u001b[39m=mu+rho*x_\u001b[39m\u001b[39m{\u001b[39m\u001b[39mt-1}+sigma*eps_\u001b[39m\u001b[39m{t}\u001b[39;00m\u001b[39m'\u001b[39m] \u001b[39m# define the law of motion of x_t\u001b[39;00m\n\u001b[0;32m <a href='vscode-notebook-cell:/c%3A/Users/MBBar/OneDrive/Escritorio/git_local/dsolve/notebooks/dsolve_examples.ipynb#X20sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m calibration \u001b[39m=\u001b[39m {\u001b[39m'\u001b[39m\u001b[39mrho\u001b[39m\u001b[39m'\u001b[39m:\u001b[39m0.8\u001b[39m,\u001b[39m'\u001b[39m\u001b[39msigma\u001b[39m\u001b[39m'\u001b[39m:\u001b[39m1\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mmu\u001b[39m\u001b[39m'\u001b[39m:\u001b[39m1\u001b[39m} \u001b[39m# give numerical values to parameters.\u001b[39;00m\n\u001b[1;32m----> <a href='vscode-notebook-cell:/c%3A/Users/MBBar/OneDrive/Escritorio/git_local/dsolve/notebooks/dsolve_examples.ipynb#X20sZmlsZQ%3D%3D?line=4'>5</a>\u001b[0m system \u001b[39m=\u001b[39m Klein(eq, x\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mx_\u001b[39m\u001b[39m{\u001b[39m\u001b[39mt-1}\u001b[39m\u001b[39m'\u001b[39m, z\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39meps_\u001b[39m\u001b[39m{t}\u001b[39;00m\u001b[39m'\u001b[39m, calibration\u001b[39m=\u001b[39mcalibration)\n\u001b[0;32m <a href='vscode-notebook-cell:/c%3A/Users/MBBar/OneDrive/Escritorio/git_local/dsolve/notebooks/dsolve_examples.ipynb#X20sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m system\u001b[39m.\u001b[39msimulate(x0 \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m, z \u001b[39m=\u001b[39m {\u001b[39m'\u001b[39m\u001b[39meps_\u001b[39m\u001b[39m{0}\u001b[39;00m\u001b[39m'\u001b[39m:\u001b[39m1\u001b[39m}, T\u001b[39m=\u001b[39m\u001b[39m12\u001b[39m)\n",
"File \u001b[1;32mc:\\users\\mbbar\\onedrive\\escritorio\\git_local\\dsolve\\src\\dsolve\\solvers.py:90\u001b[0m, in \u001b[0;36mKlein.__init__\u001b[1;34m(self, equations, x, p, z, s, indices, calibration)\u001b[0m\n\u001b[0;32m 88\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcalibrate(calibration)\n\u001b[0;32m 89\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msteady_state \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mlinalg\u001b[39m.\u001b[39minv(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39msystem_numeric[\u001b[39m'\u001b[39m\u001b[39mA\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m-\u001b[39m\u001b[39mself\u001b[39m\u001b[39m.\u001b[39msystem_numeric[\u001b[39m'\u001b[39m\u001b[39mB\u001b[39m\u001b[39m'\u001b[39m])\u001b[39m@self\u001b[39m\u001b[39m.\u001b[39msystem_numeric[\u001b[39m'\u001b[39m\u001b[39mC\u001b[39m\u001b[39m'\u001b[39m]\n\u001b[1;32m---> 90\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msystem_solution \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msolve()\n\u001b[0;32m 91\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m 92\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcalibration, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msystem_numeric, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msolution \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m, \u001b[39mNone\u001b[39;00m, \u001b[39mNone\u001b[39;00m\n",
"File \u001b[1;32mc:\\users\\mbbar\\onedrive\\escritorio\\git_local\\dsolve\\src\\dsolve\\solvers.py:219\u001b[0m, in \u001b[0;36mKlein.solve\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 216\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m'\u001b[39m\u001b[39mSystem is not calibrated.\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[0;32m 218\u001b[0m system_numeric \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msystem_numeric\n\u001b[1;32m--> 219\u001b[0m A_0, A_1, gamma \u001b[39m=\u001b[39m system_numeric[\u001b[39m'\u001b[39;49m\u001b[39mA_0\u001b[39;49m\u001b[39m'\u001b[39;49m], system_numeric[\u001b[39m'\u001b[39m\u001b[39mA_1\u001b[39m\u001b[39m'\u001b[39m], system_numeric[\u001b[39m'\u001b[39m\u001b[39mgamma\u001b[39m\u001b[39m'\u001b[39m]\n\u001b[0;32m 220\u001b[0m S, T, _, _, Q, Z \u001b[39m=\u001b[39m ordqz(A_0, A_1, output\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mcomplex\u001b[39m\u001b[39m'\u001b[39m,sort\u001b[39m=\u001b[39m\u001b[39mlambda\u001b[39;00m alpha,beta: np\u001b[39m.\u001b[39mabs(beta\u001b[39m/\u001b[39m(alpha\u001b[39m+\u001b[39m\u001b[39m1e-10\u001b[39m))\u001b[39m<\u001b[39m\u001b[39m1\u001b[39m)\n\u001b[0;32m 221\u001b[0m Q \u001b[39m=\u001b[39m Q\u001b[39m.\u001b[39mconjugate()\u001b[39m.\u001b[39mT\n",
"\u001b[1;31mKeyError\u001b[0m: 'A_0'"
]
}
],
"source": [
"eq = ['x_{t}=mu+rho*x_{t-1}+sigma*eps_{t}'] # define the law of motion of x_t\n",
"\n",
"calibration = {'rho':0.8,'sigma':1,'mu':1} # give numerical values to parameters.\n",
"\n",
"system = Klein(eq, x='x_{t-1}', z='eps_{t}', calibration=calibration)\n",
"system.simulate(x0 = 0, z = {'eps_{0}':1}, T=12)"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down
32 changes: 23 additions & 9 deletions src/dsolve/atoms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,16 @@ def is_indexed(self):
indices = self.indices
return np.any([str(i) in ['i','j','k', 'l'] for i in indices])

def reset_t(self):
'''
>>> Variable('x_{0}').reset_t()
x_{t}
'''
indices = self.indices.copy()
indices[-1]='t'
return Variable.from_elements(self.base, indices, e_t=self.e_t)


@staticmethod
def split(name:str)->tuple[Expr, Symbol, list[Expr]]:
'''
Expand All @@ -48,7 +58,7 @@ def split(name:str)->tuple[Expr, Symbol, list[Expr]]:
name = re.search('(?<=\[).*(?=\])',name).group()
base = Symbol(re.sub('_{.*?}','',name))
if re.search('(?<=_{).+?(?=})',name) is None:
raise ValueError('Variable needs to have at least one index')
raise ValueError(f'Variable {name} needs to have at least one index')
indices = re.search('(?<=_{).+?(?=})',name).group()
indices = re.split(',',indices) if ',' in indices else re.split('(?=[ijklt])',indices)
indices = [sympify(i) for i in indices if i!='']
Expand Down Expand Up @@ -138,18 +148,22 @@ def lead(self, periods:int=1):
return self.lag(-periods)

class Parameter:
'''
indices: ijkl
issues: rho_{v,i} is not allowed. rho^v_i is.
'''
def __init__(self, name):
name = normalize_string(name)
self.base = Symbol(re.sub('_{.*?}','',name))
#self.base = Symbol(re.sub('_{.*?}','',name))
self.indices = None
self.indexed = False
if re.search('(?<=_{).+?(?=})',name) is not None:
self.indexed = True
self.indexed = re.search('(?<=_{)[0-9ijkl,]+?(?=})',name) is not None
if self.indexed:
self.base = Symbol(re.sub('_{.*?}','',name))
indices = re.search('(?<=_{).+?(?=})',name).group()
indices = re.split(',',indices)
if len(indices)>1:
raise ValueError('Code cannot handle more than one index for parameters')
self.indices = [Symbol(i) for i in indices]
indices = re.split(',',indices) if ',' in indices else re.split('(?=[ijkl])',indices)
self.indices = [sympify(i) for i in indices if i!='']
else:
self.base = Symbol(name)
self.sympy = Symbol(str(self))

def subs(self, x:dict|float):
Expand Down
9 changes: 7 additions & 2 deletions src/dsolve/expressions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,20 @@
from .utils import normalize_string, normalize_dict
import re
import numpy as np
from sympy import Eq, Expr, Symbol
from sympy import Eq, Expr, Symbol, log, exp
import sympy as sym

def flatten(l):
return [item for sublist in l for item in sublist]


class DynamicExpression:
def __init__(self, expression:str):
self.elements = split(str(expression))
self.variables = {str(Variable(i)):Variable(i) for i in self.elements if is_variable(i)}
self.parameters = {str(Parameter(i)):Parameter(i) for i in self.elements if is_parameter(i)}
self.indexed = np.any([v.indexed for v in self.variables.values()])
self.indices = list(set(flatten([v.indices for v in self.variables.values()])))

@property
def sympy(self):
Expand Down Expand Up @@ -167,7 +172,7 @@ def classify_string(string):
return 'number'
elif re.search('_{[^\\\]*t.*}', string) is not None:
return 'variable'
elif string in ('+','-','/','*','(',')','='):
elif string in ('+','-','/','*','(',')','=','log','exp'):
return 'operator'
else:
return 'parameter'
Expand Down
64 changes: 64 additions & 0 deletions src/dsolve/sims.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"eq = [\n",
" 'y_t=Ey_{t+1}+\\frac{1}{\\sigma}*(v_t)',\n",
" '\\pi_t=\\beta E\\pi_{t+1}+kappa*y_t'\n",
"\n",
"]\n",
"\n",
"y= ''"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"class Sims:\n",
" def __init__(self):\n",
" pass"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.6 64-bit",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "d18f7005fc2c8fe4eb58f7d7cc5b45c2a267e69313ade3fbf0c464a72a083ce2"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
73 changes: 45 additions & 28 deletions src/dsolve/solvers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ class SystemIndices:
start:list[int]
end: list[int]

class Klein():
class Klein:

@property
def free_symbols(self):
Expand All @@ -77,6 +77,7 @@ def __init__(self,
raise ValueError(f'More unknowns ({self.n_x+self.n_p}) than equations ({self.n_eq}) ')
self.type = self.classify_system()
self.system_symbolic = self.get_matrices()

self.system_numeric = None
self.system_solution = None
if calibration is not None:
Expand Down Expand Up @@ -137,43 +138,51 @@ def read_system(self,

def expand_index(self, l, i, start, stop):
out = []

ind=list(product(*[list(range(start, end+1)) for start, end in zip(indices.start, indices.end)]))

for el in l:
if el.indexed:
out = out+ [el.subs({k:v for k,v in zip(indices.indices,i)}) for i in ind]
if sym.Symbol(i) in el.indices:
out=out+[el.subs({i:j}) for j in range(start,stop+1)]
else:
out.append(el)
return out

def expand_indices(self, l:list=None, indices:SystemIndices=None)->list[Variable]:
if indices is None or l is None:
return l
out = l.copy()
for i in range(len(indices.indices)):
out = self.expand_index(out, indices.indices[i],indices.start[i],indices.end[i])
return out

def expand_calibration(self, calibration:dict)->dict:
index = self.indices.index
start, end = self.indices.start, self.indices.end
for k,v in calibration.items():
if len(np.atleast_1d(v))>1:
calibration = calibration | {str(Parameter(k).subs({index:i})):v[i] for i in range(start, end+1)}
def expand_calibration_index(self, calibration, index):
index, start, stop = index
for k,v in calibration.items():
if Parameter(k).indices is not None and sym.Symbol(index) in Parameter(k).indices:
calibration=calibration|{str(Parameter(k).subs({index:i})):v[i-start] for i in range(start, stop+1)}
calibration.pop(k)
return calibration


def expand_calibration(self, calibration:dict)->dict:
indices = self.indices
for i in range(len(indices.indices)):
calibration=self.expand_calibration_index(calibration, (indices.indices[i], indices.start[i],indices.end[i]))
return calibration

@staticmethod
def split(string)->list[str]:
l = re.split('(?<=,)|(?=,)',string)
l = close_brackets(l)
l = [i for i in l if i!='' and i!=',']
return l


def get_matrices(self)->dict[np.ndarray]:
'''
Given the system of equations, write it in the form:
A_0y(t+1) = A_1@y(t)+gamma@z(t)
'''
A_0,_ = sym.linear_eq_to_matrix([i for i in self.equations.dynamic.symbolic], self.vars.x1+self.vars.p1)
A_1,_ = sym.linear_eq_to_matrix(_,self.vars.x+self.vars.p)
gamma, _ = sym.linear_eq_to_matrix(_, self.vars.z)
return {'A_0':A_0, 'A_1':A_1, 'gamma':-gamma}
A,_ = sym.linear_eq_to_matrix([i for i in self.equations.dynamic.symbolic], self.vars.x1+self.vars.p1)
B,_ = sym.linear_eq_to_matrix(_,self.vars.x+self.vars.p)
gamma, C = sym.linear_eq_to_matrix(_, self.vars.z)
return {'A':A, 'B':B, 'gamma':-gamma, 'C': C}

def calibrate(self, calibration: dict[float], inplace=False)->dict[np.array]:
'''
Expand All @@ -185,8 +194,12 @@ def calibrate(self, calibration: dict[float], inplace=False)->dict[np.array]:
self.equations.dynamic.calibrated = [eq.subs(calibration) for eq in self.equations.dynamic.symbolic]
self.equations.static.calibrated = [eq.subs(calibration) for eq in self.equations.static.symbolic]
self.parameters.calibration = calibration
self.system_numeric = {k: np.array(v.subs(calibration)).astype(np.float64) for k,v in self.system_symbolic.items()}

try:
self.system_numeric = {k: np.array(v.subs(calibration)).astype(np.float64) for k,v in self.system_symbolic.items()}
except:
print({str(i) for i in self.parameters()}.difference(calibration.keys()))
raise ValueError('Error')

def solve(self)->dict[np.ndarray]:
'''
Expand All @@ -200,18 +213,20 @@ def solve(self)->dict[np.ndarray]:
raise ValueError('System is not calibrated.')

system_numeric = self.system_numeric
A_0, A_1, gamma = system_numeric['A_0'], system_numeric['A_1'], system_numeric['gamma']
S, T, _, _, Q, Z = ordqz(A_0, A_1, output='complex',sort=lambda alpha,beta: np.abs(beta/(alpha+1e-10))<1)
A, B, gamma = system_numeric['A'], system_numeric['B'], system_numeric['gamma']
S, T, _, _, Q, Z = ordqz(A, B, output='complex',sort=lambda alpha,beta: np.round(np.abs(beta/np.maximum(alpha,1e-15)),6)<=1)
Q = Q.conjugate().T
#n_s = len([i for i in np.abs(np.diag(T)/np.diag(S)) if i<1]) #number of stable eigenvalues
n_s = len([np.abs(T[i,i]/S[i,i]) for i in range(S.shape[0]) if np.abs(S[i,i])>1e-6 and np.abs(T[i,i]/S[i,i])<1])
n_s = len([_ for i in range(S.shape[0]) if np.abs(S[i,i])>1e-6 and np.round(np.abs(T[i,i]/S[i,i]),6)<=1])
#print(f'System with {n_s} stable eigenvalues and {self.n_x} pre-determined variables.')

if n_s>len(self.vars.x):
raise ValueError('Multiple solutions')
print(f'Eigenvalues: {np.diag(np.abs(T/S))}')
raise ValueError(f'Multiple solutions: {n_s} stable eigenvalues for {len(self.vars.x)} pre-determined variables')

elif n_s<len(self.vars.x):
raise ValueError('No solution')
print(f'Eigenvalues: {np.diag(np.abs(T/S))}')
raise ValueError(f'No solution: {n_s} stable eigenvalues for {len(self.vars.x)} pre-determined variables')

if self.type == 'forward-looking':
return {'N': np.real(Z@(-inv(T)@Q@gamma))}
Expand Down Expand Up @@ -251,7 +266,7 @@ def normalize_z(self, z:dict, T=None)->np.array:
out = {str(k):np.zeros(T, dtype=float) for k in self.vars.z}
for iz in z:
t = Variable(iz).indices[-1]
out[f'{str(Variable(iz).base)}_{{t}}'][t]=z[str(Variable(iz))]
out[str(Variable(iz).reset_t())][t]=z[str(Variable(iz))]
else:
T = np.max([len(v) for v in z.values()])
out = {str(k):np.zeros(T,dtype=float) for k in self.vars.z}
Expand Down Expand Up @@ -337,13 +352,15 @@ def solve_static(self, d:dict[np.array])->dict[np.array]:
d[str(s.lhs)]= np.array(s_t)
return d



class MITShock:
def __init__(self, d:dict, model:Klein):
self.d = d
self.model = model

def __call__(self, var:str):
var = str(Variable(var))
return self.d[var]

def plot(self, ax, vars:str, **kwargs):

vars=[str(Variable(i)) for i in vars.split(',')]
Expand Down Expand Up @@ -373,4 +390,4 @@ def plot_expr(self, ax, eq:str, **kwargs):
if 'legend' in kwargs and kwargs['legend']:
out.append(ax.legend())

return out
return out
Loading

0 comments on commit 070a8aa

Please sign in to comment.