From 070a8aa5844b39c55bbfbd8a3bb9bb7f1e6b68b8 Mon Sep 17 00:00:00 2001 From: Marc de la Barrera i Bardalet Date: Thu, 1 Dec 2022 14:50:12 -0500 Subject: [PATCH] commit --- notebooks/dsolve_examples.ipynb | 45 +++++++- src/dsolve/atoms.py | 32 ++++-- src/dsolve/expressions.py | 9 +- src/dsolve/sims.ipynb | 64 +++++++++++ src/dsolve/solvers.py | 73 +++++++----- src/dsolve/state_space.ipynb | 193 ++++++++++++++++++++++++++++++++ src/dsolve/statespace.py | 146 ++++++++++++++++++++++++ tests/test_atoms.py | 1 + tests/test_expressions.py | 5 +- 9 files changed, 523 insertions(+), 45 deletions(-) create mode 100644 src/dsolve/sims.ipynb create mode 100644 src/dsolve/state_space.ipynb create mode 100644 src/dsolve/statespace.py diff --git a/notebooks/dsolve_examples.ipynb b/notebooks/dsolve_examples.ipynb index d19f45a..4484919 100644 --- a/notebooks/dsolve_examples.ipynb +++ b/notebooks/dsolve_examples.ipynb @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -34,7 +34,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -56,7 +56,7 @@ "E_{t}[\\pi^{w}_{t+1}]" ] }, - "execution_count": 10, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } @@ -82,7 +82,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -94,7 +94,7 @@ " 0.13421773, 0.10737418])}" ] }, - "execution_count": 11, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -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\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\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 3\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----> 5\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 6\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": {}, diff --git a/src/dsolve/atoms.py b/src/dsolve/atoms.py index 8b4c7e4..d7691c1 100644 --- a/src/dsolve/atoms.py +++ b/src/dsolve/atoms.py @@ -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]]: ''' @@ -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!=''] @@ -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): diff --git a/src/dsolve/expressions.py b/src/dsolve/expressions.py index 95bb96e..d78117f 100644 --- a/src/dsolve/expressions.py +++ b/src/dsolve/expressions.py @@ -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): @@ -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' diff --git a/src/dsolve/sims.ipynb b/src/dsolve/sims.ipynb new file mode 100644 index 0000000..5f1ab7e --- /dev/null +++ b/src/dsolve/sims.ipynb @@ -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 +} diff --git a/src/dsolve/solvers.py b/src/dsolve/solvers.py index c118593..1874a93 100644 --- a/src/dsolve/solvers.py +++ b/src/dsolve/solvers.py @@ -52,7 +52,7 @@ class SystemIndices: start:list[int] end: list[int] -class Klein(): +class Klein: @property def free_symbols(self): @@ -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: @@ -137,26 +138,35 @@ 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) @@ -164,16 +174,15 @@ def split(string)->list[str]: 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]: ''' @@ -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]: ''' @@ -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_snp.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} @@ -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(',')] @@ -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 \ No newline at end of file + return out diff --git a/src/dsolve/state_space.ipynb b/src/dsolve/state_space.ipynb new file mode 100644 index 0000000..9a4460c --- /dev/null +++ b/src/dsolve/state_space.ipynb @@ -0,0 +1,193 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "from dsolve import statespace" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "s=statespace.StateSpace(['x_{it}=0.9*x_{it-1}+z_{t}'], x='x_{it}',s='',z='z_{t}', indices={'i':(0,2)})" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "([Eq(x_{0,0}, 2.8),\n", + " Eq(x_{1,0}, 1.0),\n", + " Eq(x_{2,0}, 3.7),\n", + " Eq(x_{0,1}, 0.9*x_{0,0}),\n", + " Eq(x_{1,1}, 0.9*x_{1,0}),\n", + " Eq(x_{2,1}, 0.9*x_{2,0}),\n", + " Eq(x_{0,2}, 0.9*x_{0,1}),\n", + " Eq(x_{1,2}, 0.9*x_{1,1}),\n", + " Eq(x_{2,2}, 0.9*x_{2,1}),\n", + " Eq(x_{0,3}, 0.9*x_{0,2}),\n", + " Eq(x_{1,3}, 0.9*x_{1,2}),\n", + " Eq(x_{2,3}, 0.9*x_{2,2}),\n", + " Eq(x_{0,4}, 0.9*x_{0,3}),\n", + " Eq(x_{1,4}, 0.9*x_{1,3}),\n", + " Eq(x_{2,4}, 0.9*x_{2,3}),\n", + " Eq(x_{0,5}, 0.9*x_{0,4}),\n", + " Eq(x_{1,5}, 0.9*x_{1,4}),\n", + " Eq(x_{2,5}, 0.9*x_{2,4}),\n", + " Eq(x_{0,6}, 0.9*x_{0,5}),\n", + " Eq(x_{1,6}, 0.9*x_{1,5}),\n", + " Eq(x_{2,6}, 0.9*x_{2,5}),\n", + " Eq(x_{0,7}, 0.9*x_{0,6}),\n", + " Eq(x_{1,7}, 0.9*x_{1,6}),\n", + " Eq(x_{2,7}, 0.9*x_{2,6}),\n", + " Eq(x_{0,8}, 0.9*x_{0,7}),\n", + " Eq(x_{1,8}, 0.9*x_{1,7}),\n", + " Eq(x_{2,8}, 0.9*x_{2,7}),\n", + " Eq(x_{0,9}, 0.9*x_{0,8}),\n", + " Eq(x_{1,9}, 0.9*x_{1,8}),\n", + " Eq(x_{2,9}, 0.9*x_{2,8})],\n", + " {x_{2,2}: 2.99700000000000,\n", + " x_{2,4}: 2.42757000000000,\n", + " x_{1,6}: 0.531441000000000,\n", + " x_{1,7}: 0.478296900000000,\n", + " x_{2,7}: 1.76969853000000,\n", + " x_{1,9}: 0.387420489000000,\n", + " x_{1,0}: 1.00000000000000,\n", + " x_{2,8}: 1.59272867700000,\n", + " x_{2,1}: 3.33000000000000,\n", + " x_{0,5}: 1.65337200000000,\n", + " x_{0,0}: 2.80000000000000,\n", + " x_{0,3}: 2.04120000000000,\n", + " x_{1,5}: 0.590490000000000,\n", + " x_{0,7}: 1.33923132000000,\n", + " x_{0,8}: 1.20530818800000,\n", + " x_{0,6}: 1.48803480000000,\n", + " x_{0,9}: 1.08477736920000,\n", + " x_{1,2}: 0.810000000000000,\n", + " x_{1,3}: 0.729000000000000,\n", + " x_{1,8}: 0.430467210000000,\n", + " x_{1,1}: 0.900000000000000,\n", + " x_{2,9}: 1.43345580930000,\n", + " x_{0,4}: 1.83708000000000,\n", + " x_{0,1}: 2.52000000000000,\n", + " x_{0,2}: 2.26800000000000,\n", + " x_{2,6}: 1.96633170000000,\n", + " x_{2,3}: 2.69730000000000,\n", + " x_{1,4}: 0.656100000000000,\n", + " x_{2,0}: 3.70000000000000,\n", + " x_{2,5}: 2.18481300000000})" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "s.solve(10,z={'z_{0}':1}, x0={'x_{0,-1}':2, 'x_{1,-1}':0,'x_{2,-1}':3})" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "a=_[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "a=_[1]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{x_{2,2}: 2.99700000000000,\n", + " x_{2,4}: 2.42757000000000,\n", + " x_{1,6}: 0.531441000000000,\n", + " x_{1,7}: 0.478296900000000,\n", + " x_{2,7}: 1.76969853000000,\n", + " x_{1,9}: 0.387420489000000,\n", + " x_{1,0}: 1.00000000000000,\n", + " x_{2,8}: 1.59272867700000,\n", + " x_{2,1}: 3.33000000000000,\n", + " x_{0,5}: 1.65337200000000,\n", + " x_{0,0}: 2.80000000000000,\n", + " x_{0,3}: 2.04120000000000,\n", + " x_{1,5}: 0.590490000000000,\n", + " x_{0,7}: 1.33923132000000,\n", + " x_{0,8}: 1.20530818800000,\n", + " x_{0,6}: 1.48803480000000,\n", + " x_{0,9}: 1.08477736920000,\n", + " x_{1,2}: 0.810000000000000,\n", + " x_{1,3}: 0.729000000000000,\n", + " x_{1,8}: 0.430467210000000,\n", + " x_{1,1}: 0.900000000000000,\n", + " x_{2,9}: 1.43345580930000,\n", + " x_{0,4}: 1.83708000000000,\n", + " x_{0,1}: 2.52000000000000,\n", + " x_{0,2}: 2.26800000000000,\n", + " x_{2,6}: 1.96633170000000,\n", + " x_{2,3}: 2.69730000000000,\n", + " x_{1,4}: 0.656100000000000,\n", + " x_{2,0}: 3.70000000000000,\n", + " x_{2,5}: 2.18481300000000}" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "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 +} diff --git a/src/dsolve/statespace.py b/src/dsolve/statespace.py new file mode 100644 index 0000000..f4b2f7e --- /dev/null +++ b/src/dsolve/statespace.py @@ -0,0 +1,146 @@ + +from __future__ import annotations +import sympy as sym +from .expressions import DynamicEquation, close_brackets, DynamicExpression +from dataclasses import dataclass, field +from .utils import normalize_string, normalize_dict +import numpy as np +from .atoms import Variable, E, Parameter +import re + +@dataclass +class SystemVariables: + ''' + Class that reads and contains all information regarding system variables. + ''' + x : list[sym.Symbol] = field(default_factory=list) + z : list[sym.Symbol] = field(default_factory=list) + s : list[sym.Symbol] = field(default_factory=list) + +@dataclass +class SystemEquations(): + dynamic: SystemEquationsDynamic + static: SystemEquationsStatic = None + +@dataclass +class SystemEquationsDynamic(): + symbolic: list[sym.Eq] + calibrated: list[sym.Eq] = None + +@dataclass +class SystemEquationsStatic(): + symbolic: list[sym.Eq] + calibrated: list[sym.Eq] = None + +@dataclass +class SystemParameters: + parameters: list[sym.Symbol] + calibration: dict[float] = None + + def __call__(self): + return self.parameters + +@dataclass +class SystemIndices: + indices: list[str] + start:list[int] + end: list[int] + +class StateSpace: + def __init__(self, equations, x,s,z, indices): + indices = self.read_indices(indices) + equations = [DynamicEquation(eq) for eq in equations] + if indices is not None: + equations = self.expand_indices(equations, indices) + self.equations = equations + self.vars = self.read_variables(x,z,s,indices) + + def read_variables(self, x,z,s, indices:SystemIndices=None)->SystemVariables: + x,z,s = [self.split(i) if isinstance(i,str) else i for i in [x,z,s]] + x,z,s = [[Variable(j) for j in i] if i is not None else [] for i in [x,z,s]] + if indices is not None: + x,z,s = [self.expand_indices(i, indices) for i in [x,z,s]] + x,z,s = [[j.sympy for j in i] for i in [x,z,s]] + return SystemVariables(x,z,s) + + @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 read_equations(self, equations, s, indices): + equations = [DynamicEquation(eq) for eq in equations] + if indices is not None: + equations = self.expand_indices(equations, indices) + if s == []: + dynamic_equations =SystemEquationsDynamic([eq.sympy for eq in equations]) + static_equations = SystemEquationsStatic([]) + else: + static_equations = [eq for eq in equations if eq.lhs in s] + d = {str(eq.lhs):eq.rhs for eq in static_equations} + dynamic_equations = [eq for eq in equations if eq not in static_equations] + for i, eq in enumerate(dynamic_equations): + if len(eq.free_symbols.intersection(s))>0: + dynamic_equations[i]=eq.subs(d) + static_equations =SystemEquationsStatic([eq.sympy for eq in static_equations]) + dynamic_equations =SystemEquationsDynamic([eq.sympy for eq in dynamic_equations]) + return SystemEquations(dynamic_equations, static_equations) + + def read_indices(self, indices:dict[list[int]])->SystemIndices: + if indices is None: + return None + indices, start, end = list(indices.keys()), [indices[i][0] for i in indices], [indices[i][1] for i in indices] + indices = SystemIndices(indices, start, end) + return indices + + + def expand_indices(self, l:list=None, indices:SystemIndices=None): + 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_index(self, l, i, start, stop): + out = [] + for el in l: + 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 normalize_z(self, z:dict, T=None)->np.array: + ''' + + >>> self.normalize_z({'z_{0}':1.},T=4) + np.array([1.,0.,0.,0.]) + + >>> self.normalize_z({'z_{0}':1.,'z_2':2.},T=4) + np.array([1.,0.,2.,0.]) + ''' + z = normalize_dict(z) + out={} + if T is not None: + for t in range(T): + for iz in self.vars.z: + out=out|{str(Variable(iz).subs({'t':t})):0} + return out|z + + def solve(self, T, z, x0): + z = self.normalize_z(z,T) + system=[] + x=set() + for t in range(T): + for eq in self.equations: + eq = eq.subs({'t':t}) + eq = eq.subs(z) + eq = eq.subs(x0) + system.append(eq.sympy) + x = x.union(eq.free_symbols) + A,b = sym.linear_eq_to_matrix(system, x) + sol = A.inv()@b + return system, {k:v for k,v in zip(x,sol)} \ No newline at end of file diff --git a/tests/test_atoms.py b/tests/test_atoms.py index bc01055..58c4d44 100644 --- a/tests/test_atoms.py +++ b/tests/test_atoms.py @@ -29,6 +29,7 @@ def test_split(): def test_Parameter(): assert str(Parameter('\beta'))=='\\beta' assert str(Parameter('\beta'))==r'\beta' + assert Parameter('\rho_{v}').indexed==False assert str(Parameter('\rho_{i}'))=='\\rho_{i}' assert str(Parameter('\sigma'))=='\sigma' assert str(Parameter('\theta'))=='\\theta' diff --git a/tests/test_expressions.py b/tests/test_expressions.py index 826f73b..84d7932 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -38,4 +38,7 @@ def test_from_sympy(): def test_subs(): assert DynamicExpression('i_t-E_t[x_{t+1}]').subs({'i_t':3}) assert DynamicExpression('i_t-E_t[x_{t+1}]').subs({'i_t':3, 'E_t[x_{t+1}]':2}) - assert float(DynamicExpression('i_t-E_t[x_{t+1}]').subs({'i_t':3, 'E_t[x_{t+1}]':2}))==1. \ No newline at end of file + assert float(DynamicExpression('i_t-E_t[x_{t+1}]').subs({'i_t':3, 'E_t[x_{t+1}]':2}))==1. + +def test_indices(): + assert DynamicEquation('x_{it}=y_{it}').indices == ['i','t'] \ No newline at end of file