diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml new file mode 100644 index 00000000..b3f571d0 --- /dev/null +++ b/.github/workflows/ci.yaml @@ -0,0 +1,35 @@ +name: CI + +on: [push] + +jobs: + + build: + + runs-on: ubuntu-latest + strategy: + matrix: + python-version: [ '3.7', '3.8' ] + + name: Dolo Python Version ${{ matrix.python-version }} + steps: + + - uses: actions/checkout@v2 + + - name: Setup python + uses: actions/setup-python@v1 + with: + python-version: ${{ matrix.python-version }} + architecture: x64 + + - name: Install Poetry + uses: Gr1N/setup-poetry@v2 + + # - name: Code Quality + # run: poetry run black . --check + + - name: Install dependencies + run: poetry install + + - name: Test with pytest + run: poetry run pytest diff --git a/.github/workflows/publish-to-test-pypi.yml b/.github/workflows/publish-to-test-pypi.yml new file mode 100644 index 00000000..0626fc21 --- /dev/null +++ b/.github/workflows/publish-to-test-pypi.yml @@ -0,0 +1,41 @@ +name: Publish Python 🐍 distributions 📦 to PyPI and TestPyPI + +on: push + +jobs: + build-n-publish: + name: Build and publish Python 🐍 distributions 📦 to PyPI and TestPyPI + runs-on: ubuntu-18.04 + + + steps: + - uses: actions/checkout@master + - name: Set up Python 3.7 + uses: actions/setup-python@v1 + with: + python-version: 3.7 + - name: Install pep517 + run: >- + python -m + pip install + pep517 + --user + - name: Build a binary wheel and a source tarball + run: >- + python -m + pep517.build + --source + --binary + --out-dir dist/ + . + - name: Publish distribution 📦 to Test PyPI + uses: pypa/gh-action-pypi-publish@master + with: + password: ${{ secrets.pypitest_password }} + repository_url: https://test.pypi.org/legacy/ + + - name: Publish distribution 📦 to PyPI + if: startsWith(github.event.ref, 'refs/tags') + uses: pypa/gh-action-pypi-publish@master + with: + password: ${{ secrets.pypi_pass }} diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 64b2e352..00000000 --- a/.travis.yml +++ /dev/null @@ -1,41 +0,0 @@ -language: python - -python: - - "3.6" - -sudo: false - -#addons: -# apt: -# packages: -# - gfortran -# - liblapack-dev - -install: - # You may want to periodically update this, although the conda update - # conda line below will keep everything up-to-date. We do this - # conditionally because it saves us some downloading if the version is - # the same. - - - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; - - bash miniconda.sh -b -p $HOME/miniconda - - export PATH="$HOME/miniconda/bin:$PATH" - - hash -r - - conda config --set always_yes yes --set changeps1 no - - conda update -q conda - # Useful for debugging any issues with conda - - conda info -a - - - deps='pip nose ipython jsonschema pyyaml numba numpy scipy numexpr sympy matplotlib runipy mistune xarray' - - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION $deps - - source activate test-environment - - conda install -c conda-forge ruamel.yaml quantecon interpolation - - conda install -c albop dolang - - python setup.py install - - -script: - - - nosetests - # - runipy examples/notebooks/rbc_model.ipynb - # - runipy examples/notebooks/sudden_stops.ipynb diff --git a/BUGS b/BUGS deleted file mode 100644 index 50d6ef09..00000000 --- a/BUGS +++ /dev/null @@ -1 +0,0 @@ -- when approximating a decision rule, sigma is a numpy array of dtype object \ No newline at end of file diff --git a/COPYING b/COPYING new file mode 100644 index 00000000..f0488b0f --- /dev/null +++ b/COPYING @@ -0,0 +1,25 @@ +BSD 2-Clause License + +Copyright (c) [year], [fullname] +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/bin/dolo-julia b/bin/dolo-julia deleted file mode 100755 index 491a68a6..00000000 --- a/bin/dolo-julia +++ /dev/null @@ -1,65 +0,0 @@ -#!/usr/bin/python - -import argparse -import os - -from dolo import __version__ - -parser = argparse.ArgumentParser(description='Julia compiler') -parser.add_argument('-v','--version', action='version', version=__version__) -parser.add_argument('-r','--print_residuals', action='store_const', const=True, default=False, help='print residuals at the steady-state') -#parser.add_argument('-s','--solve', action='store_const', const=True, default=False, help='solve for the decision rule') -#parser.add_argument('-o','--order', nargs=1, type=int, default=[1], help='solution order (1,2,3)') -parser.add_argument('input', help='model file') -parser.add_argument('output',nargs='?',type=str,default=None,help='model file') - -args = parser.parse_args() - -###### - -input_file = args.input - -# note input_rad is the full path with truncated filename -[input_rad, extension] = os.path.splitext(input_file) - -if extension == '': - extension = '.yaml' -elif extension != '.yaml': - print('Unknown filetype : {}'.format(extension)) - exit(1) - -filename = input_rad + extension - -if args.output: - output_filename = args.output -else: # we should determine some good output name in case none has been specified - output_filename = input_rad + '_model.jl' - -###### - -from dolo.misc.yamlfile import yaml_import - -model = yaml_import( filename, compiler=None ) - -import re - -basename = os.path.basename(output_filename) -fname = re.compile('(.*)\.jl').match(basename).group(1) -model.name = fname - - -# check steady-state -if args.print_residuals: - from dolo.symbolic.model import print_residuals - print_residuals(model) - - -from dolo.compiler.compiler_julia import CompilerJulia -comp = CompilerJulia(model) - -txt = comp.process_output() - -###### - -with file(output_filename,'w') as f: - f.write(txt) \ No newline at end of file diff --git a/bin/dolo-matlab b/bin/dolo-matlab deleted file mode 100755 index 94c6f797..00000000 --- a/bin/dolo-matlab +++ /dev/null @@ -1,85 +0,0 @@ -#!/usr/bin/python - -import argparse -import os - -from dolo import __version__ - -parser = argparse.ArgumentParser(description='Matlab compiler') -parser.add_argument('-v','--version', action='version', version=__version__) -parser.add_argument('-d','--diff', action='store_const', const=True, default=False, help='compute symbolic derivatives') -parser.add_argument('-r','--print_residuals', action='store_const', const=True, default=False, help='print residuals at the steady-state') -parser.add_argument('-t','--model_spec', nargs=1, type=str, default=None, help='type of model') -parser.add_argument('--recipes', nargs=1, type=str, default=None, help='file containing recipes') - -parser.add_argument('input', help='model file') -parser.add_argument('output',nargs='?',type=str,default=None,help='model file') - -args = parser.parse_args() - -###### - -input_file = args.input - -# note input_rad is the full path with truncated filename -[input_rad, extension] = os.path.splitext(input_file) - -if extension == '': - extension = '.yaml' -elif extension != '.yaml': - print('Unknown filetype : {}'.format(extension)) - exit(1) - -filename = input_rad + extension - -if args.output: - output_filename = args.output -else: # we should determine some good output name in case none has been specified - output_filename = input_rad + '.m' - -###### - -from dolo import yaml_import - -model = yaml_import( filename ) - -import re - -basename = os.path.basename(output_filename) -fname = re.compile('(.*)\.m').match(basename).group(1) -model.name = fname - - -# check steady-state -if args.print_residuals: - - print(model) - -if args.recipes: - import yaml - with file(args.recipes[0]) as f: - recipes = yaml.load(f) -else: - recipes = None - -if args.diff: - diff = True -else: - diff = False - -if args.model_spec: - model_spec = args.model_spec[0] -else: - model_spec = None - -from dolo.compiler.function_compiler_matlab import compile_model_matlab - -txt = compile_model_matlab(model) -# comp = CompilerMatlab(model, recipes=recipes, model_spec=model_spec) -# -# txt = comp.process_output( diff=diff ) - -###### - -with file(output_filename,'w') as f: - f.write(txt) diff --git a/bin/dolo-recs b/bin/dolo-recs deleted file mode 100755 index f6876c2a..00000000 --- a/bin/dolo-recs +++ /dev/null @@ -1,50 +0,0 @@ -#!/usr/bin/python - -import argparse -import os -import re -from dolo import __version__ -from dolo.misc.yamlfile import yaml_import -from dolo.compiler.compiler_recs import RecsCompiler - -# Parse input arguments -parser = argparse.ArgumentParser(description='RECS compiler') -parser.add_argument('-v', '--version', action='version', version=__version__) -parser.add_argument('input', help='model file') -parser.add_argument('output', nargs='?', type=str, default=None, help='model file') - -args = parser.parse_args() - -input_file = args.input - -[input_rad, extension] = os.path.splitext(input_file) - -if extension == '': - extension = '.yaml' -elif extension != '.yaml': - print('Unknown filetype : {}'.format(extension)) - exit(1) - -filename = input_rad + extension - -# Output name -if args.output: - output_filename = args.output -else: # we should determine some good output name in case none has been specified - output_filename = input_rad + 'model' + '.m' - -# Parse yaml file -model = yaml_import(input_file, compiler=None) - -# Model name based on output file name -basename = os.path.basename(output_filename) -fname = re.compile('(.*)\.m').match(basename).group(1) -model['name'] = fname - -# Compilation for Matlab -comp = RecsCompiler(model) -txt = comp.process_output_recs() - -# Write output -with file(output_filename,'w') as f: - f.write(txt) diff --git a/bin/doloc b/bin/doloc deleted file mode 100755 index b76c7f67..00000000 --- a/bin/doloc +++ /dev/null @@ -1,213 +0,0 @@ -#! /usr/bin/python - -# To change this template, choose Tools | Templates -# and open the template in the editor. - -__author__="pablo" -__date__ ="$14 juin 2009 22:38:48$" - -import sys -import os -import getopt -#import commands - -#from dolo import * -import re -from dolo.misc.modfile import dynare_import -from dolo.compiler.compiler_dynare import DynareCompiler -from dolo.compiler.compiler_mirfac import MirFacCompiler - -if __name__ == "__main__": - - def main(argv): - if len(argv)<1: - print("Not enough arguments.") - print('') - usage() - - sys.exit(2) - - # Read options - options = {"check":False,"ramsey":False, "dynare":False, "output":False, "static-order":2, "dynamic-order":1, "mirfac": False, 'target':'recs', - 'solve': False - } - short_arg_dict = "hcdrom" - long_arg_dict = ["help","check","dynare","static-order=","dynamic-order=","output=","ramsey","mirfac","target=",'solve'] - try: - opts, args = getopt.getopt(argv, short_arg_dict, long_arg_dict) - except getopt.GetoptError: - usage() - sys.exit(2) - for opt, arg in opts: - if opt in ("-h","--help"): - usage() - sys.exit() - if opt in ("-d","--dynare"): - options["dynare"] = True - if opt in ("-c","--check"): - options["check"] = True - if opt in ("-r","--ramsey"): - options["ramsey"] = True - if opt in ("-o","--output"): - options["output"] = True - if opt in ("--dynare",): - options["dynare"] - if opt in ("--static-order",): - options["static-order"] = int(arg) - if opt in ("--dynamic-order",): - options["dynamic-order"] = int(arg) - if opt in ("--mirfac",'-m'): - options["mirfac"] = True - if opt in ("--target"): - options["target"] = arg - if opt in ("--solve"): - options["solve"] = True - print arg - - if args == []: - print("File argument missing") - sys.exit() - else: - filename = args[0] - - # determine filename type - regex_mod = re.compile("(.*)\.mod") - regex_mod_match = re.match(regex_mod,filename) - if regex_mod_match: - filetype = "mod" - filename_trunc = regex_mod_match.groups()[0] - - regex_yaml = re.compile("(.*)\.yaml") - regex_yaml_match = re.match(regex_yaml,filename) - if regex_yaml_match: - filetype = "yaml" - filename_trunc = regex_yaml_match.groups()[0] - - - current_dir = os.getcwd() - filename_trunc = current_dir + '/' + filename_trunc - - if filetype == "": - print("Unknown filetype") - sys.exit(2) - - # Import the model - if filetype == "mod": - model = dynare_import(filename) - - elif filetype == "yaml": - from dolo.misc.yamlfile import yaml_import - model = yaml_import(filename) - - - model.check_consistency(verbose=True) - - #elif filetype == "xml": - # import_xmlfile(filename_trunc,options) - - # Process the model - if options['ramsey']: - from dolo.symbolic.ramsey import RamseyModel - #dynare_model.introduce_auxiliary_variables() - model.check() - - rmodel = RamseyModel(dynare_model) - rmodel.check(verbose=True) - options['output'] = True - dynare_model=rmodel - - - if options['output']: - from dolo.compiler.compiler_dynare import DynareCompiler - comp = DynareCompiler(model) - comp.export_to_modfile() - - if options['dynare']: - from dolo.compiler.compiler_dynare import DynareCompiler - - comp = DynareCompiler(model) - - write_file(model.fname + '_dynamic.m', comp.compute_dynamic_mfile(max_order=options["dynamic-order"])) - write_file(model.fname + '_static.m', comp.compute_static_mfile(max_order=options["static-order"])) - write_file(model.fname + '.m', comp.compute_main_file() ) - - - if options['mirfac']: - from dolo.compiler.compiler_mirfac import MirFacCompiler - - mirfac_target = options['target'] - comp = MirFacCompiler(model) - - if mirfac_target == 'recs': - out_txt = comp.process_output_matlab(mirfac_target, with_parameters_values=True, with_solution=options['solve']) - else: - out_txt = comp.process_output_matlab_old(mirfac_target, with_parameters_values=True, with_solution=options['solve']) - - write_file( '{0}_{1}.m'.format(model.fname, mirfac_target) , out_txt ) - - def write_file(fname,content): - f = file(fname,'w') - f.write(content) - f.close() - - def import_modfile(filename_trunc,options): - from dolo.compiler.compiler_dynare import DynareCompiler - import time - - t0 = time.time() - '''read mod file ; write xml file''' - #from misc.interactive import dynare_import - filename = filename_trunc + ".mod" - - fname = re.compile('.*/(.*)').match(filename_trunc).group(1) - - model = dynare_import(filename) - - model.fname = fname - - model.check_consistency(verbose=True) - - comp = DynareCompiler(model) - #comp.export_infos() - - - #write_file(model.fname + '.m', comp.compute_main_file() ) - print('Modfile preprocessing finished in {0} seconds'.format(time.time() - t0)) -# if options["check"]: -# model.check_all(print_info=True,print_eq_info=True) -# if options["portfolio"]: -# from misc.portfolio import process_portfolio_model,process_portfolio_model_total -# pf_solver = process_portfolio_model(solver) -# print(pf_solver.export_to_modfile()) - - - def write_ramsey_policy(filename_trunc,options,model): - #model.export_to_modfile() - ramsey_model = model.process_ramsey_policy_model() - f = open(filename_trunc + "_ramsey.mod","w") - f.write(ramsey_model.export_to_modfile()) - f.close() - print("Ramsey model written to : " + filename_trunc + "_ramsey.mod") - - def usage(): - help_text = ''' - Usage : dolo.py [options] model.yaml|model.mod - - The file argument defining a model can be in Dynare format or in YAML format. Its extension determines its type. - - Available options are : - - -h --help print this message - -c --check model is checked for consistency - - -r --ramsey model's equations defining Ramsey's optimal policy are added (not implemented) - - -d --dynare - --static-order=order _static.m file is written at given order - --dynamic-order=order _dynamic.m file is written at given order - - -m --miranda-fackler - ''' - print(help_text) - - main(sys.argv[1:]) diff --git a/bin/dyngui b/bin/dyngui deleted file mode 100644 index f2ef17f3..00000000 --- a/bin/dyngui +++ /dev/null @@ -1,396 +0,0 @@ -#!/usr/bin/env python - -import dolo -import sys - -try: - from PyQt4.uic import loadUiType - from PyQt4 import QtCore, QtGui - from PyQt4.QtCore import QSettings - from PyQt4.QtCore import pyqtSignal as Signal - from PyQt4.QtCore import pyqtSlot as Slot - from PyQt4.QtCore import QFile - from PyQt4 import uic - engine = 'pyqt' -except Exception as e: - raise(e) - from load_ui import loadUiType - from PySide import QtCore, QtGui - from PySide.QtCore import QSettings - from PySide.QtUiTools import QUiLoader - from PySide.QtCore import QFile - engine = 'pyside' - -app = QtGui.QApplication(sys.argv) -app.setOrganizationName('Dolo') -app.setOrganizationName('Pablo') - -[EqWidgetUi,EqWidgetBase] = loadUiType("equation_widget.ui") - -class MainWindow(QtGui.QMainWindow): -# def accept(self):h -# print 'hello' -# self.eql.add_equation() - current_file = None - - def __init__(self): - QtGui.QMainWindow.__init__(self) - - # Set up the user interface from Designer. - - f = QFile("modfile_editor.ui") - f.open(QFile.ReadOnly) - - if engine == 'pyside': - loader = QUiLoader() - myWidget = loader.load(f, self) - else: - myWidget = uic.loadUi(f) - - - f.close() -# self.ui = loadUi("modfile_editor.ui") - self.ui = myWidget - -# self.ew = uic.loadUi("equation_widget.ui") - - self.ui.optionsDialog = OptionsDialog(self) - - - # Connect up the buttons. - self.connect(self.ui.pushButton, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("add_widget()")) - self.connect(self.ui.pushButton_2, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("check()")) - self.connect(self.ui.actionOpen, QtCore.SIGNAL("triggered()"), self, QtCore.SLOT("open()")) - self.connect(self.ui.actionSave, QtCore.SIGNAL("triggered()"), self, QtCore.SLOT("save_as()")) - self.connect(self.ui.actionGeneral_options, QtCore.SIGNAL("triggered()"), self.ui.optionsDialog, QtCore.SLOT("show()")) - - self.add_widget() - - self.ui.show() - - -# def restore_settings(self): -# self.settings. - - - @Slot() - def open(self): - filename = QtGui.QFileDialog.getOpenFileName() - filename = str(filename[0]) - with file(filename) as f: - txt = f.read() - try: - model = dolo.dynare_import(filename) - model.check_consistency(verbose=True) - n = len(model.equations) - q = len(self.widgets) - if n > q: - for i in range(n-q): - self.add_widget() - elif n < q: - for i in range(q-n): - self.delete_widget(0) - for n,eq in enumerate(model.equations): - tt = str(eq) - tt = tt.replace('**','^') - tt = tt.replace('==','=') - self.widgets[n].set_text(tt) - self.ui.lineEdit_var.setText( str.join(' ',[ str(v) for v in model.variables ]) ) - self.ui.lineEdit_parameters.setText( str.join(' ',[ str(v) for v in model.parameters ]) ) - self.ui.lineEdit_shocks.setText( str.join(' ',[ str(v) for v in model.shocks ]) ) - - # parameters values - parms_text = '' - for p in model.parameters: - pv = model.parameters_values[p] - parms_text += '{0} = {1};\n'.format(p,pv) - self.ui.plainTextEdit.setPlainText(parms_text) - initval_text = '' - for v in model.init_values: - vv = model.init_values[v] - initval_text += '{0} = {1};\n'.format(v,vv) - self.ui.plainTextEdit_2.setPlainText(initval_text) - - sigma = model.get_distribition().sigma - if sigma: - cov_text = '' - for i,si in enumerate(model.shocks): - for j,sj in enumerate(model.shocks): - v = model.covariances[i,j] - vt = str(v).replace('**','^') - if v==0: - pass - elif i>j: - pass - elif i0: - content += ';' - simple_mod = ''' -var {vars}; -varexo {varexo}; -parameters {parms}; -{params_definition} -model; -{equations} -end; -initval; -{initval_string} -end; - '''.format(vars = var_string,varexo = shocks_string, - parms=parameters_string,equations = content, - params_definition=params_definition, initval_string=initval_string) - - from dolo.misc.modfile import parse_dynare_text - try: - model = parse_dynare_text(simple_mod) - model.check() - info = model.info - info['name'] = model.fname - txt =''' -Model check {name}: -Number of variables : {n_variables} -Number of equations : {n_equations} -Number of shocks : {n_shocks}'''.format(**info) - except Exception as e: - msg = '\nModel is not valid' - self.ui.textEdit.setText(msg) - print(simple_mod) - print(e) - return -# try: - from dolo.numeric.perturbations import solve_decision_rule - dr = solve_decision_rule( model, order = 1 ) - txt += '\nBlanchard-Kahn conditions are met.' -# except ImportError as e: -# txt += '\nImpossible to solve the model (lapack could not be imported)' -# except Exception as e: -# txt += '\nImpossible to solve the model (yet).' -# print e - self.ui.textEdit.setText(txt) - - - def declared_symbols(self): - vars = [dolo.Variable(e) for e in str.split(str(self.ui.lineEdit_var.text()),' ' ) ] - varexos = [dolo.Shock(e) for e in str.split(str(self.ui.lineEdit_shocks.text()),' ' ) ] - parms = [dolo.Parameter(e) for e in str.split(str(self.ui.lineEdit_parameters.text()),' ') ] - context = {} - for s in vars + varexos + parms: - context[s.name] = s - return context - -import sympy - -class CustomWidget(QtGui.QWidget): - - def setup(self,n): - - self.n = n - self.ui = EqWidgetUi() - self.ui.setupUi(self) - self.ui.frame.setVisible(False) - - self.connect(self.ui.toolButton, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("specialSignal()")) - self.connect(self.ui.pushButton, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("check()")) - self.connect(self.ui.pushButton_2, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("def_parameter()")) - self.connect(self.ui.pushButton_3, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("def_variable()")) - self.connect(self.ui.pushButton_4, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("def_shock()")) - - self.locals = {} - - - @Slot() - def def_parameter(self): - t = self.undefined - tt = str( self.father.ui.lineEdit_parameters.text() ) - self.father.ui.lineEdit_parameters.setText( (tt + ' ' + t).strip(' ') ) - self.check() - - @Slot() - def def_shock(self): - t = self.undefined - tt = str( self.father.ui.lineEdit_shocks.text() ) - self.father.ui.lineEdit_shocks.setText(( tt + ' ' + t).strip(' ')) - self.check() - - @Slot() - def def_variable(self): - t = self.undefined - tt = str( self.father.ui.lineEdit_var.text() ) - self.father.ui.lineEdit_var.setText( ( tt + ' ' + t) .strip(' ') ) - self.check() - - @Slot() - def specialSignal(self): - self.emit(QtCore.SIGNAL('hi(int)'),self.n) - - @Slot() - def check(self): - txt = self.get_text() - d = self.father.declared_symbols() - known_symbols = [sympy.exp,sympy.log,sympy.sin,sympy.cos, - sympy.tan, sympy.pi,sympy.atan,sympy.sqrt] - for s in known_symbols: - d[str(s)] = s - - try: - txt = txt.replace('^','**') - if '=' in txt: - lhs,rhs = txt.split('=') - txt = lhs + '==' + rhs - eval(txt,d) - self.undefined = None - msg = 'Hiding' - self.ui.label.setText(msg) - self.ui.frame.setVisible(False) - except NameError as e: - import re - missing_name = re.compile("name \'(.*)\'").findall(str(e))[0] - msg = "Symbol ''{0}'' is not defined. Define it as :".format(missing_name) - self.undefined = missing_name - self.ui.label.setText(msg) - self.ui.frame.setVisible(True) - - - def set_number(self,n): - self.n = n - - def get_text(self): - txt = self.ui.lineEdit.text() - if txt == '': - return None - else: - return str(txt) - - def set_text(self,tt): - self.ui.lineEdit.setText(tt) - - - -[ui_class, ui_base] = loadUiType('options.ui') - - -class OptionsDialog(ui_class, ui_base): - - def __init__(self,father): - super(OptionsDialog,self).__init__() - self.setupUi(self) - self.father = father - self.connect(self.pushButton,QtCore.SIGNAL('clicked()'),self,QtCore.SLOT('open_lib_dir()')) - - ld = str(QSettings().value('lib_dir')) - if not ld: - if sys.platform == 'win32': - ld = 'c:\Windows\System32' - else: - ld = '/usr/lib/lapack' - self.set_lib_dir(ld) - - - @Slot() - def open_lib_dir(self): - lib_dir = QtGui.QFileDialog.getExistingDirectory() - print lib_dir - self.set_lib_dir( str(lib_dir) ) - - def set_lib_dir(self,lib_dir): - import ctypes, sys - try: - self.lineEdit.setText(lib_dir) - if sys.platform == 'win32': - prefix = '' - suffix = '.dll' - else: - prefix = 'lib' - suffix = '.so' - - f1 = lib_dir + '/' + prefix + 'libiomp5md' + suffix - f2 = lib_dir + '/' + prefix + 'lapack' + suffix - - - if sys.platform == 'win32': - try: - self.libiomp5md = ctypes.cdll.LoadLibrary(f1) - except Exception: - None - #self.father.lapack = - lapack = ctypes.cdll.LoadLibrary(f2) - import dolo.numeric.extern.qz - dolo.numeric.extern.qz.lapack = lapack - - QSettings().setValue("lib_dir", lib_dir) - self.label_3.setText('Loaded') - - except Exception as e: - self.label_3.setText('Not Found / Load Failed') - print e - - - - - -window = MainWindow() -sys.exit(app.exec_()) diff --git a/bin/equation_widget.ui b/bin/equation_widget.ui deleted file mode 100644 index f6c490e5..00000000 --- a/bin/equation_widget.ui +++ /dev/null @@ -1,92 +0,0 @@ - - - Form - - - - 0 - 0 - 677 - 459 - - - - Form - - - - - - - - - - - Test - - - - - - - Qt::NoFocus - - - Delete - - - - - - - - - true - - - QFrame::StyledPanel - - - QFrame::Raised - - - - - - Symbol {i} not defined - - - - - - - - - Parameter - - - - - - - Variable - - - - - - - Shock - - - - - - - - - - - - - diff --git a/bin/load_ui.py b/bin/load_ui.py deleted file mode 100644 index 3d2c9287..00000000 --- a/bin/load_ui.py +++ /dev/null @@ -1,27 +0,0 @@ -import pysideuic -import xml.etree.ElementTree as xml -from cStringIO import StringIO - -from PySide import QtGui - -def loadUiType(uiFile): - """ - Pyside lacks the "loadUiType" command, so we have to convert the ui file to py code in-memory first - and then execute it in a special frame to retrieve the form_class. - """ - parsed = xml.parse(uiFile) - widget_class = parsed.find('widget').get('class') - form_class = parsed.find('class').text - - with open(uiFile, 'r') as f: - o = StringIO() - frame = {} - - pysideuic.compileUi(f, o, indent=0) - pyc = compile(o.getvalue(), '', 'exec') - exec pyc in frame - - #Fetch the base_class and form class based on their type in the xml from designer - form_class = frame['Ui_%s'%form_class] - base_class = eval('QtGui.%s'%widget_class) - return form_class, base_class diff --git a/bin/modfile_editor.ui b/bin/modfile_editor.ui deleted file mode 100644 index 2b9cc920..00000000 --- a/bin/modfile_editor.ui +++ /dev/null @@ -1,320 +0,0 @@ - - - MainWindow - - - - 0 - 0 - 966 - 834 - - - - MainWindow - - - - - - - 0 - - - - Tab 1 - - - - - - QFrame::StyledPanel - - - QFrame::Raised - - - - - - - - Variables - - - - - - - Shocks - - - - - - - Parameters - - - - - - - Undeclared - - - - - - - - - - - - - - - - false - - - You cannot edit this field - - - - - - - - - - - - true - - - - - 0 - 0 - 376 - 462 - - - - - - - Qt::Vertical - - - - 358 - 118 - - - - - - - - - - - - - - - - - Add equation - - - - - - - Check - - - - - - - - - - Tab 2 - - - - - - - - - - 0 - 0 - 966 - 25 - - - - - File - - - - - - - - - Options - - - - - - - - - - 8 - - - - - - - - 1 - - - - - - - 2 - - - - - 0 - 0 - 98 - 96 - - - - Parameters values - - - - - - - DejaVu Sans Mono - - - - - - - - - - 0 - 0 - 98 - 96 - - - - Init values - - - - - - - DejaVu Sans Mono - - - - - - - - - - 0 - 0 - 274 - 602 - - - - Shocks definition - - - - - - - DejaVu Sans Mono - - - - - - - - - - - - - - 2 - - - - - - - - DejaVu Sans Mono - - - - - - - - - - Open - - - - - Save - - - - - Exit - - - - - General options - - - - - - diff --git a/bin/options.ui b/bin/options.ui deleted file mode 100644 index 14b50142..00000000 --- a/bin/options.ui +++ /dev/null @@ -1,115 +0,0 @@ - - - Dialog - - - - 0 - 0 - 401 - 132 - - - - Application settings - - - true - - - - - - Shared libraries - - - - - - - - Location for libs/dlls : - - - - - - - - - - Browse - - - - - - - - - - - LAPACK : - - - - - - - Not found - - - - - - - - - - - - Qt::Horizontal - - - QDialogButtonBox::Cancel|QDialogButtonBox::Ok - - - - - - - - - buttonBox - accepted() - Dialog - accept() - - - 248 - 254 - - - 157 - 274 - - - - - buttonBox - rejected() - Dialog - reject() - - - 316 - 260 - - - 286 - 274 - - - - - diff --git a/dolo-matlab.spec b/dolo-matlab.spec deleted file mode 100644 index 39cd2a52..00000000 --- a/dolo-matlab.spec +++ /dev/null @@ -1,23 +0,0 @@ -# -*- mode: python -*- -a = Analysis(['bin/dolo-matlab'], - hiddenimports=['scipy.special._ufuncs_cxx'], - hookspath=None, - runtime_hooks=None) -pyz = PYZ(a.pure) -exe = EXE(pyz, - a.scripts, - exclude_binaries=True, - name='dolo-matlab.exe', - debug=False, - strip=None, - upx=True, - console=True ) - -coll = COLLECT(exe, - a.binaries, - a.zipfiles, - a.datas, - [('recipes.yaml', 'dolo\\compiler\\recipes.yaml', 'DATA')], - strip=None, - upx=True, - name='dolo-matlab') \ No newline at end of file diff --git a/dolo/compiler/model.py b/dolo/compiler/model.py index 5bc52429..aec45f7e 100644 --- a/dolo/compiler/model.py +++ b/dolo/compiler/model.py @@ -283,13 +283,19 @@ def decode_complementarity(comp, control): class Model(SymbolicModel): - def __init__(self, data): + def __init__(self, data, check=True): self.data = data self.model_type = 'dtcc' self.__functions__ = None # self.__compile_functions__() self.set_changed() + if check: + self.calibration + self.domain + self.exogenous + self.x_bounds + self.functions def set_changed(self): self.__domain__ = None diff --git a/dolo/compiler/recipes.py b/dolo/compiler/recipes.py index 37de315c..c939c1c3 100644 --- a/dolo/compiler/recipes.py +++ b/dolo/compiler/recipes.py @@ -17,4 +17,4 @@ DATA_PATH = os.path.join(DIR_PATH, "recipes.yaml") with open(DATA_PATH) as f: - recipes = yaml.load(f) + recipes = yaml.safe_load(f) diff --git a/dolo/misc/caching.py b/dolo/misc/caching.py index 5a668f4d..1daf93ef 100644 --- a/dolo/misc/caching.py +++ b/dolo/misc/caching.py @@ -111,20 +111,20 @@ def __setitem__(self, key, value): filename = self.get_filename(key) try: with open(filename,'w') as f: - pickle.dump(value,f) + pickle.dump(value,f) except TypeError as e: raise e - + def get(self, item): import pickle - filename = self.get_filename(item) + filename = self.get_filename(item) try: with open(filename) as f: value = pickle.load(f) return value except : return None - + import collections @@ -133,11 +133,11 @@ def get(self, item): def hashable(obj): if hasattr(obj,'flatten'): # for numpy arrays return tuple( obj.flatten().tolist() ) - if isinstance(obj, collections.Hashable): + if isinstance(obj, collections.abc.Hashable): return obj - if isinstance(obj, collections.Mapping): + if isinstance(obj, collections.abc.Mapping): items = [(k,hashable(v)) for (k,v) in obj.items()] return frozenset(items) - if isinstance(obj, collections.Iterable): + if isinstance(obj, collections.abc.Iterable): return tuple([hashable(item) for item in obj]) return TypeError(type(obj)) diff --git a/dolo/tests/test_filters.py b/dolo/tests/test_filters.py deleted file mode 100644 index 8f2050b2..00000000 --- a/dolo/tests/test_filters.py +++ /dev/null @@ -1,62 +0,0 @@ -import unittest - - - -def hpfilter_dense( x, lam=1600 ): - - import numpy as np - T = x.shape[-1] - Au = np.array([ - [1, -2], - [-2, 5] - ]) - Ad = np.array([ - [5, -2], - [-2, 1] - ]) - a = np.diag( np.ones( T )*6 ) - b = np.diag( np.ones( T-1 )*(-4), 1 ) - c = np.diag( np.ones( T-2 ), 2 ) - d = np.diag( np.ones( T-1 )*(-4), -1 ) - e = np.diag( np.ones( T-2 ), -2 ) - M = a + b + c + d + e - M[0:2,0:2] = Au - M[-2:,-2:] = Ad - M *= lam - M += np.eye(T) - - if x.ndim == 1: - return np.linalg.solve(M,x) - elif x.ndim > 3: - raise Exception('HP filter is not defined for dimension >= 3.') - else: - return np.linalg.solve(M,x.T).T - - - -class FiltersTest(unittest.TestCase): - - def test_hpfilter(self): - - from dolo.numeric.filters import hp_filter - - import numpy - import numpy.random - - sig = numpy.eye( 2 ) # number of series - N = 100 # number of observations - - dumb_series = numpy.random.multivariate_normal([0,0], sig, N).T - - trend_test = hpfilter_dense(dumb_series) - - [T, cycle] = hp_filter(dumb_series) - - from numpy.testing import assert_almost_equal - - assert_almost_equal(trend_test, T) - - -if __name__ == "__main__": - - unittest.main() \ No newline at end of file diff --git a/dolo/version.py b/dolo/version.py index 20ebd867..70f8ba1d 100644 --- a/dolo/version.py +++ b/dolo/version.py @@ -1,2 +1,2 @@ -__version_info__ = ('0', '4', '9', '9') +__version_info__ = ('0', '4', '9', '10') __version__ = '.'.join(__version_info__) diff --git a/postBuild b/postBuild deleted file mode 100644 index 9914276e..00000000 --- a/postBuild +++ /dev/null @@ -1,2 +0,0 @@ -pip install -e . - diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..92aa0d81 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,31 @@ +[tool.poetry] +name = "dolo" +version = "0.4.10" +description = "Economic Modeling in Python" +authors = ["Winant Pablo "] +license="BSD-2-Clause" + +[tool.poetry.dependencies] +python = "^3.7" +interpolation = "2.1.4" +dolang = "^0.0.8" +numpy = "^1.18.3" +ipython = "^7.13.0" +pandas = "^1.0.3" +pyyaml = "^5.3.1" +matplotlib = "^3.2.1" +"ruamel.yaml" = "^0.16.10" +scipy = "^1.4.1" +quantecon = "^0.4.7" +multimethod = "^1.3" +xarray = "^0.15.1" +multipledispatch = "^0.6.0" +twine = "^3.1.1" + + +[tool.poetry.dev-dependencies] +pytest = "^5.2" + +[build-system] +requires = ["poetry>=0.12"] +build-backend = "poetry.masonry.api" diff --git a/readthedocs.yml b/readthedocs.yml deleted file mode 100644 index 52894ee6..00000000 --- a/readthedocs.yml +++ /dev/null @@ -1,6 +0,0 @@ -conda: - file: environment.yml - -python: - version: 3 - setup_py_install: true diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index c1c4cf8e..00000000 --- a/requirements.txt +++ /dev/null @@ -1,3 +0,0 @@ -dolang -multipledispatch -multimethod diff --git a/setup.py b/setup.py deleted file mode 100644 index e41a1576..00000000 --- a/setup.py +++ /dev/null @@ -1,29 +0,0 @@ -from setuptools import setup, find_packages - -# get version number -exec(open('dolo/version.py').read()) - -setup( - name="dolo", - version=__version__, - packages=find_packages(), - package_data={'dolo.compiler': ["recipes.yaml"]}, - test_suite='dolo.tests', - scripts=[ - 'bin/dolo-recs', 'bin/dolo-matlab', 'bin/dolo-julia', 'bin/dolo', - 'bin/dolo-lint' - ], - install_requires=[ - "dolang", "pyyaml", "numba", "numpy", "sympy", "scipy", "quantecon", "pandas", - "interpolation", "ruamel.yaml", "xarray", "altair", "multipledispatch", "multimethod" - ], - extras_require={ - 'interactive': ['ipython'], - 'plots': ["matplotlib"], - }, - author="Pablo Winant", - author_email="pablo.winant@gmail.com", - description='Economic modelling in Python', - license='BSD-2', - url='http://albop.github.com/dolo/', -) diff --git a/setup_windows.py b/setup_windows.py deleted file mode 100644 index b43d962e..00000000 --- a/setup_windows.py +++ /dev/null @@ -1,50 +0,0 @@ -from setuptools import setup,find_packages -import py2exe - -from dolo import __version__ - -from glob import glob - -data_files = [( - "Microsoft.VC90.CRT", - glob("C:\\Program Files\\pythonxy\\console\\Microsoft.VC90.CRT\\*.*") -)] - -excludes = ["pywin", "pywin.debugger", "pywin.debugger.dbgcon", - "pywin.dialogs", "pywin.dialogs.list", - "Tkconstants","Tkinter","tcl","zmq"] - -import sys -sys.path.append("C:\\Program Files\\pythonxy\\console\\Microsoft.VC90.CRT\\") - -setup ( - name = 'dolo', - version = __version__, - - # Declare your packages' dependencies here, for eg: - - # Fill in these to make your Egg ready for upload to - # PyPI - author = 'Pablo Winant', - author_email = '', - - url = 'www.mosphere.fr', - license = 'BSD-2', - - packages = ['dolo'], - - scripts = ['bin/dolo-matlab','bin/dolo-recs'], - - data_files = data_files, - - options = { - "py2exe" : { - #"includes" : ["sip"], - "excludes" : excludes, - "dist_dir" : "../windist/", - "bundle_files" : 3 - } - }, - console = ['bin/dolo-matlab','bin/dolo-recs'], - -) diff --git a/trash/dolo/algos/dtcscc/__init__.py b/trash/dolo/algos/dtcscc/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/trash/dolo/algos/dtcscc/accuracy.py b/trash/dolo/algos/dtcscc/accuracy.py deleted file mode 100644 index 0426fdd7..00000000 --- a/trash/dolo/algos/dtcscc/accuracy.py +++ /dev/null @@ -1,243 +0,0 @@ -from __future__ import division - -import numpy as np - -from dolo.algos.dtcscc.simulations import simulate -from dolo.numeric.discretization.quadrature import gauss_hermite_nodes -from dolo.numeric.misc import mlinspace - - -class EulerErrors(dict): - - @property - def max_errors(self): - return self['max_errors'] - - @property - def mean_errors(self): - return self['ergodic'] - - @property - def time_weighted(self): - return self['time_weighted'] - - def __repr__(self): - measures = ['max_errors', 'ergodic'] - measures += ['time_weighted'] if 'time_weighted' in self else [] - s = 'Euler Errors:\n' - for m in measures: - s += '- {:15}: {}\n'.format(m, self[m]) - return s - - -class DenHaanErrors(dict): - - @property - def max_errors(self): - return self['max_errors'] - - @property - def mean_errors(self): - return self['mean_errors'] - - def __repr__(self): - - measures = ['max_errors', 'mean_errors'] - s = 'Den Haan errors:\n' - for m in measures: - s += '- {:15}: {}\n'.format(m, self[m]) - return s - - -def omega(model, dr, n_exp=10000, grid={}, bounds=None, - n_draws=100, seed=0, horizon=50, s0=None, - solve_expectations=False, time_discount=None): - - assert(model.model_type == 'dtcscc') - - f = model.functions['arbitrage'] - g = model.functions['transition'] - - distrib = model.get_distribution() - sigma = distrib.sigma - - parms = model.calibration['parameters'] - - mean = np.zeros(sigma.shape[0]) - - np.random.seed(seed) - epsilons = np.random.multivariate_normal(mean, sigma, n_draws) - weights = np.ones(epsilons.shape[0])/n_draws - - approx = model.get_endo_grid(**grid) - a = approx.a - b = approx.b - orders = approx.orders - bounds = np.row_stack([a, b]) - - domain = RectangularDomain(a, b, orders) - - grid = domain.grid - - n_s = len(model.symbols['states']) - - errors = test_residuals(grid, dr, f, g, parms, epsilons, weights) - errors = abs(errors) - - if s0 is None: - s0 = model.calibration['states'] - - simul = simulate(model, dr, s0, n_exp=n_exp, horizon=horizon+1, - discard=True, solve_expectations=solve_expectations, return_array=True) - - s_simul = simul[:, :, :n_s] - - densities = [domain.compute_density(s_simul[t, :, :]) - for t in range(horizon)] - ergo_dens = densities[-1] - - max_error = np.max(errors, axis=0) # maximum errors - ergo_error = np.dot(ergo_dens, errors) # weighted by erg. distr. - - d = dict( - errors=errors, - densities=densities, - bounds=bounds, - max_errors=max_error, - ergodic=ergo_error, - domain=domain - ) - - if time_discount is not None: - beta = time_discount - time_weighted_errors = max_error*0 - for i in range(horizon): - err = np.dot(densities[i], errors) - time_weighted_errors += beta**i * err - time_weighted_errors /= (1-beta**(horizon-1))/(1-beta) - d['time_weighted'] = time_weighted_errors - - return EulerErrors(d) - - -def denhaanerrors(model, dr, s0=None, horizon=100, n_sims=10, seed=0, - integration_orders=None): - - assert(model.model_type == 'dtcscc') - - n_x = len(model.symbols['controls']) - n_s = len(model.symbols['states']) - - distrib = model.get_distribution() - sigma = distrib.sigma - mean = sigma[0, :]*0 - - if integration_orders is None: - integration_orders = [5]*len(mean) - [nodes, weights] = gauss_hermite_nodes(integration_orders, sigma) - - if s0 is None: - s0 = model.calibration['states'] - - # standard simulation - simul = simulate(model, dr, s0, horizon=horizon, n_exp=n_sims, seed=seed, - solve_expectations=False, return_array=True) - simul_se = simulate(model, dr, s0, horizon=horizon, n_exp=n_sims, - seed=seed, solve_expectations=True, nodes=nodes, - weights=weights, return_array=True) - - x_simul = simul[:, n_s:n_s+n_x, :] - x_simul_se = simul_se[:, n_s:n_s+n_x, :] - - diff = abs(x_simul_se - x_simul) - error_1 = (diff).max(axis=0).mean(axis=1) - error_2 = (diff).mean(axis=0).mean(axis=1) - - d = dict( - max_errors=error_1, - mean_errors=error_2, - horizon=horizon, - n_sims=n_sims - ) - - return DenHaanErrors(d) - - -class RectangularDomain: - - def __init__(self, a, b, orders): - self.d = len(a) - self.a = a - self.b = b - self.bounds = np.row_stack([a, b]) - self.orders = np.array(orders, dtype=int) - nodes = [np.linspace(a[i], b[i], orders[i]) - for i in range(len(orders))] - - self.nodes = nodes - self.grid = mlinspace(a, b, orders) - - def find_cell(self, x): - """ - @param x: Nxd array - @return: Nxd array with line i containing the indices of cell - containing x[i, :] - """ - - inf = self.a - sup = self.b - N = x.shape[0] - indices = np.zeros((N, self.d), dtype=int) - for i in range(self.d): - xi = (x[:, i] - inf[i])/(sup[i]-inf[i]) - ni = np.floor(xi*self.orders[i]) - ni = np.minimum(np.maximum(ni, 0), self.orders[i]-1) - indices[:, i] = ni - - return np.ravel_multi_index(indices.T, self.orders) - - def compute_density(self, x): - - import time - t1 = time.time() - - cell_indices = self.find_cell(x) - t2 = time.time() - - keep = np.isfinite(cell_indices) - cell_linear_indices = cell_indices[keep] - - npoints = cell_indices.shape[0] - - counts = np.bincount(cell_linear_indices, - minlength=self.orders.prod()) - - dens = counts/npoints - - t3 = time.time() - - return dens - - -# TODO: this logic is repeated at least here and in time_iteration.py -def test_residuals(s, dr, f, g, parms, epsilons, weights): - - n_draws = epsilons.shape[0] - - x = dr(s) - - [N, n_x] = x.shape - ss = np.tile(s, (n_draws, 1)) - xx = np.tile(x, (n_draws, 1)) - ee = np.repeat(epsilons, N, axis=0) - - ssnext = g(ss, xx, ee, parms) - xxnext = dr(ssnext) - - val = f(ss, xx, ee, ssnext, xxnext, parms) - - res = np.zeros((N, n_x)) - for i in range(n_draws): - res += weights[i] * val[N*i:N*(i+1), :] - - return res diff --git a/trash/dolo/algos/dtcscc/gssa.py b/trash/dolo/algos/dtcscc/gssa.py deleted file mode 100644 index 3513ac6c..00000000 --- a/trash/dolo/algos/dtcscc/gssa.py +++ /dev/null @@ -1,190 +0,0 @@ -import time -import warnings - -import numpy as np -from numba import jit -from scipy.linalg import lstsq -from dolo.algos.dtcscc.perturbations import approximate_controls -from dolo.algos.dtcscc.simulations import simulate -from interpolation.complete_poly import (CompletePolynomial, - _complete_poly_impl, _complete_poly_impl_vec, complete_polynomial, - n_complete) - - -def gssa(model, maxit=100, tol=1e-8, dr0=None, verbose=False, - n_sim=10000, deg=3, damp=0.1, seed=42): - """ - Sketch of algorithm: - - 0. Choose levels for the initial states and the simulation length (n_sim) - 1. Obtain an initial decision rule -- here using first order perturbation - 2. Draw a sequence of innovations epsilon - 3. Iterate on the following steps: - - Use the epsilons, initial states, and proposed decision rule to - simulate model forward. Will leave us with time series of states and - controls - - Evaluate expectations using quadrature - - Use direct response to get alternative proposal for controls - - Regress updated controls on the simulated states to get proposal - coefficients. New coefficients are convex combination of previous - coefficients and proposal coefficients. Weights controlled by damp, - where damp is the weight on the old coefficients. This should be - fairly low to increase chances of convergence. - - Check difference between the simulated series of controls and the - direct response version of controls - - """ - # verify input arguments - if deg < 0 or deg > 5: - raise ValueError("deg must be in [1, 5]") - - if damp < 0 or damp > 1: - raise ValueError("damp must be in [0, 1]") - - t1 = time.time() - - # extract model functions and parameters - g = model.__original_functions__['transition'] - g_gu = model.__original_gufunctions__['transition'] - h_gu = model.__original_gufunctions__['expectation'] - d_gu = model.__original_gufunctions__['direct_response'] - p = model.calibration['parameters'] - n_s = len(model.symbols["states"]) - n_x = len(model.symbols["controls"]) - n_z = len(model.symbols["expectations"]) - n_eps = len(model.symbols["shocks"]) - s0 = model.calibration["states"] - x0 = model.calibration["controls"] - - # construct initial decision rule if not supplied - if dr0 is None: - drp = approximate_controls(model) - else: - drp = dr0 - - # set up quadrature weights and nodes - distrib = model.get_distribution() - nodes, weights = distrib.discretize() - - # draw sequence of innovations - np.random.seed(seed) - distrib = model.get_distribution() - sigma = distrib.sigma - epsilon = np.random.multivariate_normal(np.zeros(n_eps), sigma, n_sim) - - # simulate initial decision rule and do initial regression for coefs - init_sim = simulate(model, drp, horizon=n_sim, return_array=True, - forcing_shocks=epsilon) - s_sim = init_sim[:, 0, 0:n_s] - x_sim = init_sim[:, 0, n_s:n_s + n_x] - Phi_sim = complete_polynomial(s_sim.T, deg).T - coefs = np.ascontiguousarray(lstsq(Phi_sim, x_sim)[0]) - - # NOTE: the ascontiguousarray above was needed for numba to compile the - # `np.dot` in the simulation function in no python mode. Appearantly - # the array returned from lstsq is not C-contiguous - - # allocate for simulated series of expectations and next period states - z_sim = np.empty((n_sim, n_z)) - S = np.empty_like(s_sim) - X = np.empty_like(x_sim) - H = np.empty_like(z_sim) - new_x = np.empty_like(x_sim) - - # set initial states and controls - s_sim[0, :] = s0 - x_sim[0, :] = x0 - - Phi_t = np.empty(n_complete(n_s, deg)) # buffer array for simulation - - # create jitted function that will simulate states and controls, using - # the epsilon shocks from above (define here as closure over all data - # above). - @jit(nopython=True) - def simulate_states_controls(s, x, Phi_t, coefs): - for t in range(1, n_sim): - g(s[t - 1, :], x[t - 1, :], epsilon[t, :], p, s[t, :]) - - # fill Phi_t with new complete poly version of s[t, :] - _complete_poly_impl_vec(s[t, :], deg, Phi_t) - - # do inner product to get new controls - x[t, :] = Phi_t @coefs - - it = 0 - err = 10.0 - err_0 = 10 - - if verbose: - headline = '|{0:^4} | {1:10} | {2:8} | {3:8} |' - headline = headline.format('N', ' Error', 'Gain', 'Time') - stars = '-' * len(headline) - print(stars) - print(headline) - print(stars) - - # format string for within loop - fmt_str = '|{0:4} | {1:10.3e} | {2:8.3f} | {3:8.3f} |' - - while err > tol and it <= maxit: - t_start = time.time() - - # simulate with new coefficients - simulate_states_controls(s_sim, x_sim, Phi_t, coefs) - - # update expectations of z - # update_expectations(s_sim, x_sim, z_sim, Phi_sim) - z_sim[:, :] = 0.0 - for i in range(weights.shape[0]): - e = nodes[i, :] # extract nodes - # evaluate future states at each node (stores in S) - g_gu(s_sim, x_sim, e, p, S) - - # evaluate future controls at each future state - _complete_poly_impl(S.T, deg, Phi_sim.T) - np.dot(Phi_sim, coefs, out=X) - - # compute expectation (stores in H) - h_gu(S, X, p, H) - z_sim += weights[i] * H - - # get controls on the simulated points from direct_resposne - # (stores in new_x) - d_gu(s_sim, z_sim, p, new_x) - - # update basis matrix and do regression of new_x on s_sim to get - # updated coefficients - _complete_poly_impl(s_sim.T, deg, Phi_sim.T) - new_coefs = np.ascontiguousarray(lstsq(Phi_sim, new_x)[0]) - - # check whether they differ from the preceding guess - err = (abs(new_x - x_sim).max()) - - # update the series of controls and coefficients - x_sim[:, :] = new_x - coefs = (1 - damp) * new_coefs + damp * coefs - - if verbose: - # update error and print if `verbose` - err_SA = err / err_0 - err_0 = err - t_finish = time.time() - elapsed = t_finish - t_start - if verbose: - print(fmt_str.format(it, err, err_SA, elapsed)) - - it += 1 - - if it == maxit: - warnings.warn(UserWarning("Maximum number of iterations reached")) - - # compute final fime and do final printout if `verbose` - t2 = time.time() - if verbose: - print(stars) - print('Elapsed: {} seconds.'.format(t2 - t1)) - print(stars) - - cp = CompletePolynomial(deg, len(s0)) - cp.fit_values(s_sim, x_sim) - return cp diff --git a/trash/dolo/algos/dtcscc/nonlinearsystem.py b/trash/dolo/algos/dtcscc/nonlinearsystem.py deleted file mode 100644 index 6bfe48ef..00000000 --- a/trash/dolo/algos/dtcscc/nonlinearsystem.py +++ /dev/null @@ -1,187 +0,0 @@ -import numpy -import scipy.sparse -import time -from numba import jit -from dolo.numeric.serial_operations import serial_multiplication as smult -from dolo.algos.dtcscc.perturbations import approximate_controls -from dolo.algos.dtcscc.time_iteration import create_interpolator - - -def nonlinear_system(model, dr0=None, maxit=10, tol=1e-8, grid={}, distribution={}, verbose=True): - - ''' - Finds a global solution for ``model`` by solving one large system of equations - using a simple newton algorithm. - - Parameters - ---------- - model: NumericModel - "dtcscc" model to be solved - verbose: boolean - if True, display iterations - dr0: decision rule - initial guess for the decision rule - maxit: int - maximum number of iterationsd - tol: tolerance criterium for successive approximations - grid: grid options - distribution: distribution options - - Returns - ------- - decision rule : - approximated solution - ''' - - - if verbose: - headline = '|{0:^4} | {1:10} | {2:8} |' - headline = headline.format('N', ' Error', 'Time') - stars = '-'*len(headline) - print(stars) - print(headline) - print(stars) - # format string for within loop - fmt_str = '|{0:4} | {1:10.3e} | {2:8.3f} |' - - f = model.functions['arbitrage'] - g = model.functions['transition'] - p = model.calibration['parameters'] - - distrib = model.get_distribution(**distribution) - nodes, weights = distrib.discretize() - - approx = model.get_endo_grid(**grid) - ms = create_interpolator(approx, approx.interpolation) - - grid = ms.grid - - if dr0 is None: - dr = approximate_controls(model) - else: - dr = dr0 - - ms.set_values(dr(grid)) - - x = dr(grid) - x0 = x.copy() - - it = 0 - err = 10 - - - a0 = x0.copy().reshape((x0.shape[0]*x0.shape[1],)) - a = a0.copy() - - while err > tol and it < maxit: - - it += 1 - t1 = time.time() - - r, da = residuals(f, g, grid, a.reshape(x0.shape), ms, nodes, weights, p, diff=True)[:2] - r = r.flatten() - - err = abs(r).max() - - t2 = time.time() - - if verbose: - print(fmt_str.format(it, err, t2-t1)) - - if err > tol: - a -= scipy.sparse.linalg.spsolve(da, r) - - if verbose: - print(stars) - - return ms - - -@jit -def serial_to_full(m): - # m is a N * n_x * n_x array - # it is converted into a (N*n_x)*(N*n_x) sparse array M - # such that nonzero elements are M[n,i,n,j] = m[n,i,j] - - N, n_x, n_xx = m.shape - assert(n_x == n_xx) - val = numpy.zeros((N*n_x*N)) - ind_i = numpy.zeros((N*n_x*N)) - ind_j = numpy.zeros((N*n_x*N)) - t = 0 - for n in range(N): - for i in range(n_x): - for j in range(n_x): - val[t] = m[n, i, j] - ind_i[t] = n_x*n + i - ind_j[t] = n_x*n + j - t += 1 - - mat = scipy.sparse.coo_matrix((val, (ind_i, ind_j)), shape=(N*n_x, N*n_x)) - mmat = mat.tocsr() - mmat.eliminate_zeros() - return mmat - - -@jit -def diag_to_full(m): - # m is a N * n_x * N array - # it is converted into a (N*n_x)*(N*n_x) sparse array M - # such that nonzero elements are M[p,i,q,i] = m[p,i,q] - N, n_x, NN = m.shape - assert(N == NN) - val = numpy.zeros((N*n_x*N)) - ind_i = numpy.zeros((N*n_x*N)) - ind_j = numpy.zeros((N*n_x*N)) - - t = 0 - for n in range(N): - for i in range(n_x): - for nn in range(N): - val[t] = m[n, i, nn] - ind_i[t] = n_x*n + i - ind_j[t] = n_x*nn + i - t += 1 - - mat = scipy.sparse.coo_matrix((val, (ind_i, ind_j)), shape=(N*n_x, N*n_x)) - mmat = mat.tocsr() - mmat.eliminate_zeros() - return mmat - - - -def residuals(f, g, s, x, dr, nodes, weights, p, diff=True): - - N, n_x = x.shape - - dr.set_values(x) - - output = numpy.zeros((N, n_x)) - - if not diff: - for i in range(nodes.shape[0]): - E = nodes[i, :][None, :].repeat(N, axis=0) - S = g(s, x, E, p) - X = dr.interpolate(S) - t = f(s, x, E, S, X, p) - output += weights[i]*t - return output - - if diff: - - output_x = scipy.sparse.csr_matrix((N*n_x, N*n_x)) - - for i in range(nodes.shape[0]): - E = nodes[i, :][None, :].repeat(N, axis=0) - S, S_s, S_x, S_E = g(s, x, E, p, diff=True) - X, X_S, X_x = dr.interpolate(S, deriv=True, deriv_X=True) - R, R_s, R_x, R_E, R_S, R_X = f(s, x, E, S, X, p, diff=True) - - output += weights[i]*R - - t1 = weights[i]*(R_x + smult(R_S, S_x) + smult(R_X, smult(X_S, S_x))) # N.n_x.n_x - t1 = serial_to_full(t1) - t2 = weights[i] * serial_to_full(R_X) @ diag_to_full(X_x) - output_x += t1 + t2 - - return output, output_x diff --git a/trash/dolo/algos/dtcscc/parameterized_expectations.py b/trash/dolo/algos/dtcscc/parameterized_expectations.py deleted file mode 100644 index 2b9ad35e..00000000 --- a/trash/dolo/algos/dtcscc/parameterized_expectations.py +++ /dev/null @@ -1,345 +0,0 @@ -import time -import numpy as np -from dolo.algos.dtcscc.perturbations import approximate_controls -from dolo.numeric.optimize.ncpsolve import ncpsolve -from dolo.numeric.optimize.newton import (SerialDifferentiableFunction, - serial_newton) -from dolo.numeric.interpolation import create_interpolator - - -def parameterized_expectations(model, verbose=False, dr0=None, - pert_order=1, with_complementarities=True, - grid={}, distribution={}, - maxit=100, tol=1e-8, inner_maxit=100, - direct=False): - - ''' - Find global solution for ``model`` via parameterized expectations. - Controls must be expressed as a direct function of equilibrium objects. - Algorithm iterates over the expectations function in the arbitrage equation. - - Parameters: - ---------- - - model : NumericModel - ``dtcscc`` model to be solved - - verbose : boolean - if True, display iterations - - dr0 : decision rule - initial guess for the decision rule - - pert_order : {1} - if no initial guess is supplied, the perturbation solution at order - ``pert_order`` is used as initial guess - - grid : grid options - - distribution : distribution options - - maxit : maximum number of iterations - - tol : tolerance criterium for successive approximations - - inner_maxit : maximum number of iteration for inner solver - - direct : if True, solve with direct method. If false, solve indirectly - - Returns - ------- - - decision rule : - approximated solution - - ''' - - def vprint(t): - if verbose: - print(t) - - g = model.functions['transition'] - h = model.functions['expectation'] - f = model.functions['arbitrage_exp'] # f(s, x, z, p, out) - parms = model.calibration['parameters'] - if direct is True: - d = model.functions['direct_response'] - - approx = model.get_endo_grid(**grid) - grid = approx.grid - interp_type = approx.interpolation - dr = create_interpolator(approx, interp_type) # Interp for control - expect = create_interpolator(approx, interp_type) # Interp for expectation - - distrib = model.get_distribution(**distribution) - nodes, weights = distrib.discretize() - - N = grid.shape[0] - - if dr0 is None: - if pert_order == 1: - dr0 = approximate_controls(model) - - if pert_order > 1: - raise Exception("Perturbation order > 1 not supported (yet).") - - # Use initial decision rule to find initial expectation function - x_0 = dr0(grid) - x_0 = x_0.real # just in case ... - - z_0 = np.zeros((N, len(model.symbols['expectations']))) - z_new = np.zeros((N, len(model.symbols['expectations']))) - - xxnext = np.zeros((x_0.shape[0], x_0.shape[1], weights.shape[0])) - for i in range(weights.shape[0]): - e = nodes[i, :] - ssnext = g(grid, x_0, e, parms) - xxnext[:, :, i] = dr0(ssnext) - z_0 += weights[i]*h(ssnext, xxnext[:, :, i], parms) - - t1 = time.time() - it = 0 - err = 10 - err_0 = 10 - - verbit = True if verbose == 'full' else False - - if with_complementarities is True: - lbfun = model.functions['controls_lb'] - ubfun = model.functions['controls_ub'] - lb = lbfun(grid, parms) - ub = ubfun(grid, parms) - else: - lb = None - ub = None - - if verbose: - headline = '|{0:^4} | {1:10} | {2:8} | {3:8} |' - headline = headline.format('N', ' Error', 'Gain', 'Time') - stars = '-'*len(headline) - print(stars) - print(headline) - print(stars) - - # format string for within loop - fmt_str = '|{0:4} | {1:10.3e} | {2:8.3f} | {3:8.3f} |' - - while err > tol and it <= maxit: - - it += 1 - t_start = time.time() - - # update interpolation object with current guess for expectation - expect.set_values(z_0) - - xxnext[...] = 0 - if direct is True: - # Use control as direct function of arbitrage equation - xx = d(grid, expect(grid), parms) - for i in range(weights.shape[0]): - e = nodes[i, :] - ssnext = g(grid, xx, e, parms) - xxnext[:, :, i] = d(ssnext, expect(ssnext), parms) - - if with_complementarities is True: - xx = np.minimum(xx, ub) - xx = np.maximum(xx, lb) - for i in range(weights.shape[0]): - xxnext[:, :, i] = np.minimum(xxnext[:, :, i], ub) - xxnext[:, :, i] = np.maximum(xxnext[:, :, i], lb) - else: - # Find control by solving arbitrage equation - def fun(x): return f(grid, x, expect(grid), parms) - sdfun = SerialDifferentiableFunction(fun) - - if with_complementarities is True: - [xx, nit] = ncpsolve(sdfun, lb, ub, x_0, verbose=verbit, - maxit=inner_maxit) - dr.set_values(xx) # Update decision rule object - for i in range(weights.shape[0]): - e = nodes[i, :] - ssnext = g(grid, xx, e, parms) - xxnext[:, :, i] = dr(ssnext) - - # Make sure x_t+1 is within bounds - xxnext[:, :, i] = np.minimum(xxnext[:, :, i], ub) - xxnext[:, :, i] = np.maximum(xxnext[:, :, i], lb) - else: - [xx, nit] = serial_newton(sdfun, x_0, verbose=verbit) - dr.set_values(xx) # Update decision rule object - for i in range(weights.shape[0]): - e = nodes[i, :] - ssnext = g(grid, xx, e, parms) - xxnext[:, :, i] = dr(ssnext) - - # Compute the new expectation function - z_new[...] = 0 - for i in range(weights.shape[0]): - e = nodes[i, :] - ssnext = g(grid, xx, e, parms) - z_new += weights[i]*h(ssnext, xxnext[:, :, i], parms) - - # update error - err = (abs(z_new - z_0).max()) - - # Update guess for expectations function and the decision rule - z_0 = z_new.copy() - x_0 = xx - - # print error infomation if `verbose` - err_SA = err/err_0 - err_0 = err - t_finish = time.time() - elapsed = t_finish - t_start - if verbose: - print(fmt_str.format(it, err, err_SA, elapsed)) - - if it == maxit: - import warnings - warnings.warn(UserWarning("Maximum number of iterations reached")) - - # compute final fime and do final printout if `verbose` - t2 = time.time() - if verbose: - print(stars) - print('Elapsed: {} seconds.'.format(t2 - t1)) - print(stars) - - # Interpolation for the decision rule - dr.set_values(x_0) - - return dr - -import time -import numpy as np -from dolo.algos.dtcscc.perturbations import approximate_controls -from dolo.numeric.interpolation import create_interpolator - - -def parameterized_expectations_direct(model, verbose=False, dr0=None, - pert_order=1, grid={}, distribution={}, - maxit=100, tol=1e-8): - ''' - Finds a global solution for ``model`` using parameterized expectations - function. Requires the model to be written with controls as a direct - function of the model objects. - - The algorithm iterates on the expectations function in the arbitrage - equation. It follows the discussion in section 9.9 of Miranda and - Fackler (2002). - - Parameters - ---------- - model : NumericModel - "dtcscc" model to be solved - verbose : boolean - if True, display iterations - dr0 : decision rule - initial guess for the decision rule - pert_order : {1} - if no initial guess is supplied, the perturbation solution at order - ``pert_order`` is used as initial guess - grid: grid options - distribution: distribution options - maxit: maximum number of iterations - tol: tolerance criterium for successive approximations - - Returns - ------- - decision rule : - approximated solution - ''' - - t1 = time.time() - g = model.functions['transition'] - d = model.functions['direct_response'] - h = model.functions['expectation'] - parms = model.calibration['parameters'] - - if dr0 is None: - if pert_order == 1: - dr0 = approximate_controls(model) - - if pert_order > 1: - raise Exception("Perturbation order > 1 not supported (yet).") - - approx = model.get_endo_grid(**grid) - grid = approx.grid - interp_type = approx.interpolation - dr = create_interpolator(approx, interp_type) - expect = create_interpolator(approx, interp_type) - - distrib = model.get_distribution(**distribution) - nodes, weights = distrib.discretize() - - N = grid.shape[0] - z = np.zeros((N, len(model.symbols['expectations']))) - - x_0 = dr0(grid) - x_0 = x_0.real # just in case ... - h_0 = h(grid, x_0, parms) - - it = 0 - err = 10 - err_0 = 10 - - if verbose: - headline = '|{0:^4} | {1:10} | {2:8} | {3:8} |' - headline = headline.format('N', ' Error', 'Gain', 'Time') - stars = '-'*len(headline) - print(stars) - print(headline) - print(stars) - - # format string for within loop - fmt_str = '|{0:4} | {1:10.3e} | {2:8.3f} | {3:8.3f} |' - - while err > tol and it <= maxit: - - it += 1 - t_start = time.time() - - # dr.set_values(x_0) - expect.set_values(h_0) - - z[...] = 0 - for i in range(weights.shape[0]): - e = nodes[i, :] - S = g(grid, x_0, e, parms) - # evaluate expectation over the future state - z += weights[i]*expect(S) - - # TODO: check that control is admissible - new_x = d(grid, z, parms) - new_h = h(grid, new_x, parms) - - # update error - err = (abs(new_h - h_0).max()) - - # Update guess for decision rule and expectations function - x_0 = new_x.copy() - h_0 = new_h.copy() - - # print error information if `verbose` - err_SA = err/err_0 - err_0 = err - t_finish = time.time() - elapsed = t_finish - t_start - if verbose: - print(fmt_str.format(it, err, err_SA, elapsed)) - - if it == maxit: - import warnings - warnings.warn(UserWarning("Maximum number of iterations reached")) - - # compute final fime and do final printout if `verbose` - t2 = time.time() - if verbose: - print(stars) - print('Elapsed: {} seconds.'.format(t2 - t1)) - print(stars) - - # Interpolation for the decision rule - dr.set_values(x_0) - - return dr diff --git a/trash/dolo/algos/dtcscc/perfect_foresight.py b/trash/dolo/algos/dtcscc/perfect_foresight.py deleted file mode 100644 index edba7b71..00000000 --- a/trash/dolo/algos/dtcscc/perfect_foresight.py +++ /dev/null @@ -1,389 +0,0 @@ -import numpy as np -import pandas as pd -from numpy import array, atleast_2d, linspace, zeros -from scipy.optimize import root - -from dolo.algos.dtcscc.steady_state import find_deterministic_equilibrium -from dolo.numeric.optimize.ncpsolve import ncpsolve -from dolo.numeric.optimize.newton import newton -from dolo.numeric.serial_operations import serial_multiplication as smult - -def _shocks_to_epsilons(model, shocks, T): - """ - Helper function to support input argument `shocks` being one of many - different data types. Will always return a `T, n_e` matrix. - """ - n_e = len(model.calibration['shocks']) - - # if we have a DataFrame, convert it to a dict and rely on the method below - if isinstance(shocks, pd.DataFrame): - shocks = {k: shocks[k].tolist() for k in shocks.columns} - - # handle case where shocks might be a dict. Be careful to handle case where - # value arrays are not the same length - if isinstance(shocks, dict): - epsilons = np.zeros((T+1, n_e)) - for (i, k) in enumerate(model.symbols["shocks"]): - if k in shocks: - this_shock = shocks[k] - epsilons[:len(this_shock)-1, i] = this_shock[1:] - epsilons[(len(this_shock)-1):, i] = this_shock[-1] - else: - # otherwise set to value in calibration - epsilons[:, i] = model.calibration["shocks"][i] - - return epsilons - - # read from calibration if not given - if shocks is None: - shocks = model.calibration["shocks"] - - # now we just assume that shocks is array-like and try using the output of - # np.asarray(shocks) - shocks = np.asarray(shocks) - shocks = shocks.reshape((-1, n_e)) - - # until last period, exogenous shock takes its last value - epsilons = np.zeros((T+1, n_e)) - epsilons[:(shocks.shape[0]-1), :] = shocks[1:, :] - epsilons[(shocks.shape[0]-1):, :] = shocks[-1:, :] - - return epsilons - - -def deterministic_solve(model, shocks=None, start_states=None, T=100, - ignore_constraints=False, maxit=100, - initial_guess=None, verbose=False, tol=1e-6): - """ - Computes a perfect foresight simulation using a stacked-time algorithm. - - The initial state is specified either by providing a series of exogenous - shocks and assuming the model is initially in equilibrium with the first - value of the shock, or by specifying an initial value for the states. - - Parameters - ---------- - model : NumericModel - "dtcscc" model to be solved - shocks : array-like, dict, or pandas.DataFrame - A specification of the shocks to the model. Can be any of the - following (note by "declaration order" below we mean the order - of `model.symbols["shocks"]`): - - - A 1d numpy array-like specifying a time series for a single - shock, or all shocks stacked into a single array. - - A 2d numpy array where each column specifies the time series - for one of the shocks in declaration order. This must be an - `N` by number of shocks 2d array. - - A dict where keys are strings found in - `model.symbols["shocks"]` and values are a time series of - values for that shock. For model shocks that do not appear in - this dict, the shock is set to the calibrated value. Note - that this interface is the most flexible as it allows the user - to pass values for only a subset of the model shocks and it - allows the passed time series to be of different lengths. - - A DataFrame where columns map shock names into time series. - The same assumptions and behavior that are used in the dict - case apply here - - If nothing is given here, `shocks` is set equal to the - calibrated values found in `model.calibration["shocks"]` for - all periods. - - If the length of any time-series in shocks is less than `T` - (see below) it is assumed that that particular shock will - remain at the final given value for the duration of the - simulaiton. - start_states : ndarray or dict - a vector with the value of initial states, or a calibration - dictionary with the initial values of states and controls - T : int - horizon for the perfect foresight simulation - maxit : int - maximum number of iteration for the nonlinear solver - verbose : boolean - if True, the solver displays iterations - tol : float - stopping criterium for the nonlinear solver - ignore_constraints : bool - if True, complementarity constraints are ignored. - - Returns - ------- - pandas dataframe - a dataframe with T+1 observations of the model variables along the - simulation (states, controls, auxiliaries). The first observation is - the steady-state corresponding to the first value of the shocks. The - simulation should return to a steady-state corresponding to the last - value of the exogenous shocks. - - """ - - # TODO: - - # if model.model_spec == 'fga': - # from dolo.compiler.converter import GModel_fg_from_fga - # model = GModel_fg_from_fga(model) - - # definitions - n_s = len(model.calibration['states']) - n_x = len(model.calibration['controls']) - - epsilons = _shocks_to_epsilons(model, shocks, T) - - # final initial and final steady-states consistent with exogenous shocks - if start_states is None: - start_states = model.calibration - - if isinstance(start_states, dict): - # at least that part is clear - start_equilibrium = start_states - start_s = start_equilibrium['states'] - start_x = start_equilibrium['controls'] - final_s = start_equilibrium['states'] - final_x = start_equilibrium['controls'] - elif isinstance(start_states, np.ndarray): - start_s = start_states - start_x = model.calibration['controls'] - final_s = model.calibration['states'] - final_x = model.calibration['controls'] - - # if start_constraints: - # # we ignore start_constraints - # start_dict.update(start_constraints) - # final_equilibrium = start_constraints.copy() - # else: - # final_eqm = find_deterministic_equilibrium(model, - # constraints=final_dict) - # final_s = final_eqm['states'] - # final_x = final_eqm['controls'] - # - # start_s = start_states - # start_x = final_x - - # TODO: for start_x, it should be possible to use first order guess - - final = np.concatenate([final_s, final_x]) - start = np.concatenate([start_s, start_x]) - - if verbose is True: - print("Initial states : {}".format(start_s)) - print("Final controls : {}".format(final_x)) - - p = model.calibration['parameters'] - - if initial_guess is None: - initial_guess = np.row_stack([start*(1-l) + final*l - for l in linspace(0.0, 1.0, T+1)]) - - else: - if isinstance(initial_guess, pd.DataFrame): - initial_guess = np.array(initial_guess).T.copy() - initial_guess = initial_guess[:, :n_s+n_x] - initial_guess[0, :n_s] = start_s - initial_guess[-1, n_s:] = final_x - - sh = initial_guess.shape - - if model.x_bounds and not ignore_constraints: - initial_states = initial_guess[:, :n_s] - [lb, ub] = [u(initial_states, p) for u in model.x_bounds] - lower_bound = initial_guess*0 - np.inf - lower_bound[:, n_s:] = lb - upper_bound = initial_guess*0 + np.inf - upper_bound[:, n_s:] = ub - test1 = max(lb.max(axis=0) - lb.min(axis=0)) - test2 = max(ub.max(axis=0) - ub.min(axis=0)) - if test1 > 0.00000001 or test2 > 0.00000001: - msg = "Not implemented: perfect foresight solution requires that " - msg += "controls have constant bounds." - raise Exception(msg) - else: - ignore_constraints = True - lower_bound = None - upper_bound = None - - nn = sh[0]*sh[1] - - def fobj(vec): - o = det_residual(model, vec.reshape(sh), start_s, final_x, epsilons)[0] - return o.ravel() - - if not ignore_constraints: - def ff(vec): - return det_residual(model, vec.reshape(sh), start_s, final_x, - epsilons, jactype='sparse') - - x0 = initial_guess.ravel() - sol, nit = ncpsolve(ff, lower_bound.ravel(), upper_bound.ravel(), - initial_guess.ravel(), verbose=verbose, - maxit=maxit, tol=tol, jactype='sparse') - - sol = sol.reshape(sh) - - else: - - def ff(vec): - return det_residual(model, vec.reshape(sh), start_s, final_x, - epsilons, diff=False).ravel() - - x0 = initial_guess.ravel() - sol = root(ff, x0, jac=False) - res = ff(sol.x) - sol = sol.x.reshape(sh) - - if 'auxiliary' in model.functions: - colnames = (model.symbols['states'] + model.symbols['controls'] + - model.symbols['auxiliaries']) - # compute auxiliaries - y = model.functions['auxiliary'](sol[:, :n_s], sol[:, n_s:], p) - sol = np.column_stack([sol, y]) - else: - colnames = model.symbols['states'] + model.symbols['controls'] - - sol = np.column_stack([sol, epsilons]) - colnames = colnames + model.symbols['shocks'] - - ts = pd.DataFrame(sol, columns=colnames) - return ts - - -def det_residual(model, guess, start, final, shocks, diff=True, - jactype='sparse'): - ''' - Computes the residuals, the derivatives of the stacked-time system. - :param model: an fga model - :param guess: the guess for the simulated values. An `(n_s.n_x) x N` array, - where n_s is the number of states, - n_x the number of controls, and `N` the length of the simulation. - :param start: initial boundary condition (initial value of the states) - :param final: final boundary condition (last value of the controls) - :param shocks: values for the exogenous shocks - :param diff: if True, the derivatives are computes - :return: a list with two elements: - - an `(n_s.n_x) x N` array with the residuals of the system - - a `(n_s.n_x) x N x (n_s.n_x) x N` array representing the jacobian of - the system - ''' - - # TODO: compute a sparse derivative and ensure the solvers can deal with it - - n_s = len(model.symbols['states']) - n_x = len(model.symbols['controls']) - - n_e = len(model.symbols['shocks']) - N = guess.shape[0] - - p = model.calibration['parameters'] - - f = model.functions['arbitrage'] - g = model.functions['transition'] - - vec = guess[:-1, :] - vec_f = guess[1:, :] - - s = vec[:, :n_s] - x = vec[:, n_s:] - S = vec_f[:, :n_s] - X = vec_f[:, n_s:] - - e = shocks[:-1, :] - E = shocks[1:, :] - - if diff: - SS, SS_s, SS_x, SS_e = g(s, x, e, p, diff=True) - R, R_s, R_x, R_e, R_S, R_X = f(s, x, E, S, X, p, diff=True) - else: - SS = g(s, x, e, p) - R = f(s, x, E, S, X, p) - - res_s = SS - S - res_x = R - - res = np.zeros((N, n_s+n_x)) - - res[1:, :n_s] = res_s - res[:-1, n_s:] = res_x - - res[0, :n_s] = - (guess[0, :n_s] - start) - res[-1, n_s:] = - (guess[-1, n_s:] - guess[-2, n_s:]) - - if not diff: - return res - else: - - sparse_jac = False - if not sparse_jac: - - # we compute the derivative matrix - res_s_s = SS_s - res_s_x = SS_x - - # next block is probably very inefficient - jac = np.zeros((N, n_s+n_x, N, n_s+n_x)) - for i in range(N-1): - jac[i, n_s:, i, :n_s] = R_s[i, :, :] - jac[i, n_s:, i, n_s:] = R_x[i, :, :] - jac[i, n_s:, i+1, :n_s] = R_S[i, :, :] - jac[i, n_s:, i+1, n_s:] = R_X[i, :, :] - jac[i+1, :n_s, i, :n_s] = SS_s[i, :, :] - jac[i+1, :n_s, i, n_s:] = SS_x[i, :, :] - jac[i+1, :n_s, i+1, :n_s] = -np.eye(n_s) - # jac[i,n_s:,i,:n_s] = R_s[i,:,:] - # jac[i,n_s:,i,n_s:] = R_x[i,:,:] - # jac[i+1,n_s:,i,:n_s] = R_S[i,:,:] - # jac[i+1,n_s:,i,n_s:] = R_X[i,:,:] - # jac[i,:n_s,i+1,:n_s] = SS_s[i,:,:] - # jac[i,:n_s,i+1,n_s:] = SS_x[i,:,:] - # jac[i+1,:n_s,i+1,:n_s] = -np.eye(n_s) - jac[0, :n_s, 0, :n_s] = - np.eye(n_s) - jac[-1, n_s:, -1, n_s:] = - np.eye(n_x) - jac[-1, n_s:, -2, n_s:] = + np.eye(n_x) - nn = jac.shape[0]*jac.shape[1] - res = res.ravel() - jac = jac.reshape((nn, nn)) - - if jactype == 'sparse': - from scipy.sparse import csc_matrix, csr_matrix - jac = csc_matrix(jac) - # scipy bug ? I don't get the same with csr - - return [res, jac] - - -if __name__ == '__main__': - - # this example computes the response of the rbc economy to a series of - # expected productivity shocks. investment is bounded by an exogenous value - # 0.2, so that investment is constrained in the first periods - - # TODO: propose a meaningful economic example - - from dolo import yaml_import - - m = yaml_import("../../../examples/models/compat/Figv4_1191.yaml") - T = 100 - g_list = [0.2]*10+[0.4] - - # first try using a list - sol1 = deterministic_solve(m, shocks=g_list) - - # then try using a 1d array - sol2 = deterministic_solve(m, shocks=np.asarray(g_list)) - - # then try using a 2d array - g_shock = np.array(g_list)[:, None] - sol3 = deterministic_solve(m, shocks=g_shock) - - # now try using a dict - sol4 = deterministic_solve(m, shocks={"g": g_list}) - - # now try using a DataFrame - sol5 = deterministic_solve(m, shocks=pd.DataFrame({"g": g_list})) - - # check that they are all the same - for s in [sol2, sol3, sol4, sol5]: - assert max(abs(sol1-s).max()) == 0.0 - - m2 = yaml_import("../../../examples/models/compat/rmt3_ch11.yaml") - sol = deterministic_solve(m, shocks={"g": [0.2]*10+[0.4]}, T=T) diff --git a/trash/dolo/algos/dtcscc/perturbations.py b/trash/dolo/algos/dtcscc/perturbations.py deleted file mode 100644 index cd08d643..00000000 --- a/trash/dolo/algos/dtcscc/perturbations.py +++ /dev/null @@ -1,185 +0,0 @@ -import numpy as np -from numpy import column_stack, dot, eye, row_stack, zeros -from numpy.linalg import solve - -from dolo.algos.dtcscc.steady_state import find_deterministic_equilibrium -from dolo.numeric.extern.qz import qzordered -from dolo.numeric.taylor_expansion import TaylorExpansion as CDR - - -class GeneralizedEigenvaluesError(Exception): - - def __init__(self, A=None, B=None, eigval=None, cutoff=None, n_states=None): - - self.A = A - self.B = B - self.eigval = eigval - self.eigval_sorted = sorted(eigval) - self.cutoff = cutoff - self.n_states = n_states - -class GeneralizedEigenvaluesDefinition(GeneralizedEigenvaluesError): - - def __str__(self): - - return "Generalized Eigenvalues imply a 0/0 division. Undefined solution." - - -class GeneralizedEigenvaluesSelectionError(GeneralizedEigenvaluesError): - - def __str__(self): - - return "Impossible to select the {} bigger eigenvalues in a unique way.".format(self.n_states) - - - -def approximate_controls(model, verbose=False, steady_state=None, eigmax=1.0-1e-6, - solve_steady_state=False, order=1): - """Compute first order approximation of optimal controls - - Parameters: - ----------- - - model: NumericModel - Model to be solved - - verbose: boolean - If True: displays number of contracting eigenvalues - - steady_state: ndarray - Use supplied steady-state value to compute the approximation. - The routine doesn't check whether it is really a solution or not. - - solve_steady_state: boolean - Use nonlinear solver to find the steady-state - - orders: {1} - Approximation order. (Currently, only first order is supported). - - Returns: - -------- - - TaylorExpansion: - Decision Rule for the optimal controls around the steady-state. - - """ - - if order > 1: - raise Exception("Not implemented.") - - f = model.functions['arbitrage'] - g = model.functions['transition'] - - if steady_state is not None: - calib = steady_state - else: - calib = model.calibration - - if solve_steady_state: - calib = find_deterministic_equilibrium(model) - - p = calib['parameters'] - s = calib['states'] - x = calib['controls'] - e = calib['shocks'] - - distrib = model.get_distribution() - sigma = distrib.sigma - - l = g(s, x, e, p, diff=True) - [junk, g_s, g_x, g_e] = l[:4] # [el[0,...] for el in l[:4]] - - l = f(s, x, e, s, x, p, diff=True) - [res, f_s, f_x, f_e, f_S, f_X] = l # [el[0,...] for el in l[:6]] - - n_s = g_s.shape[0] # number of controls - n_x = g_x.shape[1] # number of states - n_e = g_e.shape[1] - n_v = n_s + n_x - - A = row_stack([ - column_stack([eye(n_s), zeros((n_s, n_x))]), - column_stack([-f_S , -f_X ]) - ]) - - B = row_stack([ - column_stack([g_s, g_x]), - column_stack([f_s, f_x]) - ]) - - - [S, T, Q, Z, eigval] = qzordered(A, B, 1.0-1e-8) - - Q = Q.real # is it really necessary ? - Z = Z.real - - diag_S = np.diag(S) - diag_T = np.diag(T) - - tol_geneigvals = 1e-10 - - try: - ok = sum((abs(diag_S) < tol_geneigvals) * - (abs(diag_T) < tol_geneigvals)) == 0 - assert(ok) - except Exception as e: - raise GeneralizedEigenvaluesError(diag_S=diag_S, diag_T=diag_T) - - if max(eigval[:n_s]) >= 1 and min(eigval[n_s:]) < 1: - # BK conditions are met - pass - else: - eigval_s = sorted(eigval, reverse=True) - ev_a = eigval_s[n_s-1] - ev_b = eigval_s[n_s] - cutoff = (ev_a - ev_b)/2 - if not ev_a > ev_b: - raise GeneralizedEigenvaluesSelectionError( - A=A, B=B, eigval=eigval, cutoff=cutoff, - diag_S=diag_S, diag_T=diag_T, n_states=n_s - ) - import warnings - if cutoff > 1: - warnings.warn("Solution is not convergent.") - else: - warnings.warn("There are multiple convergent solutions. The one with the smaller eigenvalues was selected.") - [S, T, Q, Z, eigval] = qzordered(A, B, cutoff) - - Z11 = Z[:n_s, :n_s] - # Z12 = Z[:n_s, n_s:] - Z21 = Z[n_s:, :n_s] - # Z22 = Z[n_s:, n_s:] - # S11 = S[:n_s, :n_s] - # T11 = T[:n_s, :n_s] - - # first order solution - # P = (solve(S11.T, Z11.T).T @ solve(Z11.T, T11.T).T) - C = solve(Z11.T, Z21.T).T - Q = g_e - - s = s.ravel() - x = x.ravel() - - A = g_s + g_x @ C - B = g_e - - dr = CDR([s, x, C]) - dr.A = A - dr.B = B - dr.sigma = sigma - - return dr - -if __name__ == '__main__': - - from dolo import * - from os.path import abspath, dirname, join - - dolo_dir = dirname(dirname(dirname(dirname(abspath(__file__))))) - model = yaml_import(join(dolo_dir, "examples", "models", "rbc.yaml")) - model.set_calibration(dumb=1) - - from dolo.algos.dtcscc.perturbations import approximate_controls - - dr = approximate_controls(model) - print(dr) diff --git a/trash/dolo/algos/dtcscc/perturbations_higher_order.py b/trash/dolo/algos/dtcscc/perturbations_higher_order.py deleted file mode 100644 index e69de29b..00000000 diff --git a/trash/dolo/algos/dtcscc/time_iteration_2.py b/trash/dolo/algos/dtcscc/time_iteration_2.py deleted file mode 100644 index 6c8604dc..00000000 --- a/trash/dolo/algos/dtcscc/time_iteration_2.py +++ /dev/null @@ -1,312 +0,0 @@ -import time -import numpy as np -from dolo.algos.dtcscc.perturbations import approximate_controls -from dolo.numeric.optimize.ncpsolve import ncpsolve -from dolo.numeric.optimize.newton import (SerialDifferentiableFunction, - serial_newton) -from dolo.numeric.interpolation import create_interpolator - - -def parameterized_expectations(model, verbose=False, dr0=None, - pert_order=1, with_complementarities=True, - grid={}, distribution={}, - maxit=100, tol=1e-8, inner_maxit=10, - direct=False): - - ''' - Find global solution for ``model`` via parameterized expectations. - Controls must be expressed as a direct function of equilibrium objects. - Algorithm iterates over the expectations function in the arbitrage equation. - - Parameters: - ---------- - - model : NumericModel - ``dtcscc`` model to be solved - - verbose : boolean - if True, display iterations - - dr0 : decision rule - initial guess for the decision rule - - pert_order : {1} - if no initial guess is supplied, the perturbation solution at order - ``pert_order`` is used as initial guess - - grid : grid options - - distribution : distribution options - - maxit : maximum number of iterations - - tol : tolerance criterium for successive approximations - - inner_maxit : maximum number of iteration for inner solver - - direct : if True, solve with direct method. If false, solve indirectly - - Returns - ------- - - decision rule : - approximated solution - - ''' - - t1 = time.time() - - g = model.functions['transition'] - h = model.functions['expectation'] - d = model.functions['direct_response'] - f = model.functions['arbitrage_exp'] # f(s, x, z, p, out) - parms = model.calibration['parameters'] - - if dr0 is None: - if pert_order == 1: - dr0 = approximate_controls(model) - - if pert_order > 1: - raise Exception("Perturbation order > 1 not supported (yet).") - - approx = model.get_endo_grid(**grid) - grid = approx.grid - interp_type = approx.interpolation - dr = create_interpolator(approx, interp_type) - expect = create_interpolator(approx, interp_type) - - distrib = model.get_distribution(**distribution) - nodes, weights = distrib.discretize() - - N = grid.shape[0] - z = np.zeros((N, len(model.symbols['expectations']))) - - x_0 = dr0(grid) - x_0 = x_0.real # just in case ... - h_0 = h(grid, x_0, parms) - - it = 0 - err = 10 - err_0 = 10 - - verbit = True if verbose == 'full' else False - - if with_complementarities is True: - lbfun = model.functions['controls_lb'] - ubfun = model.functions['controls_ub'] - lb = lbfun(grid, parms) - ub = ubfun(grid, parms) - else: - lb = None - ub = None - - if verbose: - headline = '|{0:^4} | {1:10} | {2:8} | {3:8} |' - headline = headline.format('N', ' Error', 'Gain', 'Time') - stars = '-'*len(headline) - print(stars) - print(headline) - print(stars) - - # format string for within loop - fmt_str = '|{0:4} | {1:10.3e} | {2:8.3f} | {3:8.3f} |' - - while err > tol and it <= maxit: - - it += 1 - t_start = time.time() - - # dr.set_values(x_0) - expect.set_values(h_0) - - # evaluate expectation over the future state - z[...] = 0 - for i in range(weights.shape[0]): - e = nodes[i, :] - S = g(grid, x_0, e, parms) - z += weights[i]*expect(S) - - if direct is True: - # Use control as direct function of arbitrage equation - new_x = d(grid, z, parms) - if with_complementarities is True: - new_x = np.minimum(new_x, ub) - new_x = np.maximum(new_x, lb) - else: - # Find control by solving arbitrage equation - def fun(x): return f(grid, x, z, parms) - sdfun = SerialDifferentiableFunction(fun) - - if with_complementarities is True: - [new_x, nit] = ncpsolve(sdfun, lb, ub, x_0, verbose=verbit, - maxit=inner_maxit) - else: - [new_x, nit] = serial_newton(sdfun, x_0, verbose=verbit) - - new_h = h(grid, new_x, parms) - - # update error - err = (abs(new_h - h_0).max()) - - # Update guess for decision rule and expectations function - x_0 = new_x - h_0 = new_h - - # print error infomation if `verbose` - err_SA = err/err_0 - err_0 = err - t_finish = time.time() - elapsed = t_finish - t_start - if verbose: - print(fmt_str.format(it, err, err_SA, elapsed)) - - if it == maxit: - import warnings - warnings.warn(UserWarning("Maximum number of iterations reached")) - - # compute final fime and do final printout if `verbose` - t2 = time.time() - if verbose: - print(stars) - print('Elapsed: {} seconds.'.format(t2 - t1)) - print(stars) - - # Interpolation for the decision rule - dr.set_values(x_0) - - return dr - - -import time -import numpy as np -from dolo.algos.dtcscc.perturbations import approximate_controls -from dolo.numeric.interpolation import create_interpolator - - -def parameterized_expectations_direct(model, verbose=False, dr0=None, - pert_order=1, grid={}, distribution={}, - maxit=100, tol=1e-8): - ''' - Finds a global solution for ``model`` using parameterized expectations - function. Requires the model to be written with controls as a direct - function of the model objects. - - The algorithm iterates on the expectations function in the arbitrage - equation. It follows the discussion in section 9.9 of Miranda and - Fackler (2002). - - Parameters - ---------- - model : NumericModel - "dtcscc" model to be solved - verbose : boolean - if True, display iterations - dr0 : decision rule - initial guess for the decision rule - pert_order : {1} - if no initial guess is supplied, the perturbation solution at order - ``pert_order`` is used as initial guess - grid: grid options - distribution: distribution options - maxit: maximum number of iterations - tol: tolerance criterium for successive approximations - - Returns - ------- - decision rule : - approximated solution - ''' - - t1 = time.time() - g = model.functions['transition'] - d = model.functions['direct_response'] - h = model.functions['expectation'] - parms = model.calibration['parameters'] - - if dr0 is None: - if pert_order == 1: - dr0 = approximate_controls(model) - - if pert_order > 1: - raise Exception("Perturbation order > 1 not supported (yet).") - - approx = model.get_endo_grid(**grid) - grid = approx.grid - interp_type = approx.interpolation - dr = create_interpolator(approx, interp_type) - expect = create_interpolator(approx, interp_type) - - distrib = model.get_distribution(**distribution) - nodes, weights = distrib.discretize() - - N = grid.shape[0] - z = np.zeros((N, len(model.symbols['expectations']))) - - x_0 = dr0(grid) - x_0 = x_0.real # just in case ... - h_0 = h(grid, x_0, parms) - - it = 0 - err = 10 - err_0 = 10 - - if verbose: - headline = '|{0:^4} | {1:10} | {2:8} | {3:8} |' - headline = headline.format('N', ' Error', 'Gain', 'Time') - stars = '-'*len(headline) - print(stars) - print(headline) - print(stars) - - # format string for within loop - fmt_str = '|{0:4} | {1:10.3e} | {2:8.3f} | {3:8.3f} |' - - while err > tol and it <= maxit: - - it += 1 - t_start = time.time() - - # dr.set_values(x_0) - expect.set_values(h_0) - - z[...] = 0 - for i in range(weights.shape[0]): - e = nodes[i, :] - S = g(grid, x_0, e, parms) - # evaluate expectation over the future state - z += weights[i]*expect(S) - - # TODO: check that control is admissible - new_x = d(grid, z, parms) - new_h = h(grid, new_x, parms) - - # update error - err = (abs(new_h - h_0).max()) - - # Update guess for decision rule and expectations function - x_0 = new_x - h_0 = new_h - - # print error information if `verbose` - err_SA = err/err_0 - err_0 = err - t_finish = time.time() - elapsed = t_finish - t_start - if verbose: - print(fmt_str.format(it, err, err_SA, elapsed)) - - if it == maxit: - import warnings - warnings.warn(UserWarning("Maximum number of iterations reached")) - - # compute final fime and do final printout if `verbose` - t2 = time.time() - if verbose: - print(stars) - print('Elapsed: {} seconds.'.format(t2 - t1)) - print(stars) - - # Interpolation for the decision rule - dr.set_values(x_0) - - return dr diff --git a/trash/dolo/algos/dynare/__init__.py b/trash/dolo/algos/dynare/__init__.py deleted file mode 100644 index ab7cfb08..00000000 --- a/trash/dolo/algos/dynare/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -# from perturbations import solve_decision_rule -# from simulations import impulse_response_function, stoch_simul diff --git a/trash/dolo/algos/dynare/perturbations.py b/trash/dolo/algos/dynare/perturbations.py deleted file mode 100644 index d76842cf..00000000 --- a/trash/dolo/algos/dynare/perturbations.py +++ /dev/null @@ -1,808 +0,0 @@ -from dolo.numeric.matrix_equations import second_order_solver, solve_sylvester -from dolo.numeric.tensor import sdot,mdot - - - -import numpy as np - -TOL = 1E-10 -TOL_RES = 1E-8 - -def solve_decision_rule(model,order=1,method='default',mlab=None,steady_state = None, solve_ss = False): - - - Sigma = model.get_distribution().sigma - y = model.calibration['variables'] - x = model.calibration['shocks'] - parms = model.calibration['parameters'] - - - if steady_state != None: - y0 = steady_state - else: - y0 = y - - # Build the approximation point to take the derivatives - yy = np.concatenate([y,y,y,x]) - parms = np.array(parms) - - if order==1: - f_dynamic = model.functions['f_dynamic'] - else: - f_dynamic = model.__get_compiled_functions__(order=order)['f_dynamic'] - - derivatives = f_dynamic(yy,parms, order=order) - - derivatives_ss = None - - if (method == 'sigma2'): - n_s = len(model.shocks) - n_v = len(model.variables) - nt = n_v*3 +n_s - nn = nt+1 - - if order == 2: - [f_0,f_1,f_2] = derivatives - ft_1 = f_1[:,:nt] - ft_2 = f_2[:,:nt,:nt] - f_ss = f_2[:,nt,nt] # derivatives with sigma^2 - derivatives = [f_0,ft_1,ft_2] - derivatives_ss = [f_ss] - - - elif order == 3: - [f_0,f_1,f_2,f_3] = derivatives - ft_1 = f_1[:,:nt] - ft_2 = f_2[:,:nt,:nt] - ft_3 = f_3[:,:nt,:nt,:nt] - f_ss = f_2[:,nt,nt] # derivatives with sigma^2 - f_1ss = f_3[:,:nt,nt,nt] - derivatives = [f_0,ft_1,ft_2,ft_3] - derivatives_ss = [f_ss,f_1ss] - - if method in ('sigma2','default'): - d = perturb_solver(derivatives, Sigma, max_order=order, derivatives_ss=derivatives_ss,mlab=mlab) - - elif method == 'full': - n_s = len(model.shocks) - n_v = len(model.variables) - n_p = len(model.parameters) - d = new_solver_with_p(derivatives,(n_v,n_s,n_p),max_order=order) - return d - - d['ys'] = np.array( y ) - d['Sigma'] = Sigma - - from dolo.numeric.decision_rules import DynareDecisionRule as DDR - ddr = DDR( d , model ) - return ddr - # from dolo.numeric.decision_rules import TaylorExpansion - # return TaylorExpansion(d) - - -def perturb_solver(derivatives, Sigma, max_order=2, derivatives_ss=None, mlab=None): - - - if max_order == 1: - [f_0,f_1] = derivatives - elif max_order == 2: - [f_0,f_1,f_2] = derivatives - elif max_order == 3: - [f_0,f_1,f_2,f_3] = derivatives - else: - raise Exception('Perturbations not implemented at order {0}'.format(max_order)) - derivs = derivatives - - f = derivs - n = f[0].shape[0] # number of variables - s = f[1].shape[1] - 3*n - [n_v,n_s] = [n,s] - - - f1_A = f[1][:,:n] - f1_B = f[1][:,n:(2*n)] - f1_C = f[1][:,(2*n):(3*n)] - f1_D = f[1][:,(3*n):] - - ## first order - [ev,g_x] = second_order_solver(f1_A,f1_B,f1_C) - - res = np.dot(f1_A,np.dot(g_x,g_x)) + np.dot(f1_B,g_x) + f1_C - - mm = np.dot(f1_A, g_x) + f1_B - - g_u = - np.linalg.solve( mm , f1_D ) - - if max_order == 1: - d = {'ev':ev, 'g_a': g_x, 'g_e': g_u} - return d - - # we need it for higher order - V_a = np.concatenate( [np.dot(g_x,g_x),g_x,np.eye(n_v),np.zeros((s,n))] ) - V_e = np.concatenate( [np.dot(g_x,g_u),g_u,np.zeros((n_v,n_s)),np.eye(n_s)] ) - - # Translation - - f_1 = f[1] - f_2 = f[2] - f_d = f1_A - f_a = f1_B - g_a = g_x - g_e = g_u - n_a = n_v - n_e = n_s - - # Once for all ! - A = f_a + sdot(f_d,g_a) - B = f_d - C = g_a - A_inv = np.linalg.inv(A) - - - ################## - # Automatic code # - ################## - - #----------Computing order 2 - - order = 2 - - #--- Computing derivatives ('a', 'a') - - K_aa = + mdot(f_2,[V_a,V_a]) - L_aa = np.zeros( (n_v, n_v, n_v) ) - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aa + sdot(f_d,L_aa) - if mlab==None: - g_aa = solve_sylvester(A,B,C,D) - else: - n_d = D.ndim - 1 - n_v = C.shape[1] - CC = np.kron(C,C) - DD = D.reshape( n_v, n_v**n_d ) - [err,E] = mlab.gensylv(2,A,B,C,DD,nout=2) - g_aa = - E.reshape((n_v,n_v,n_v)) # check that - is correct - - - if order < max_order: - Y = L_aa + mdot(g_a,[g_aa]) + mdot(g_aa,[g_a,g_a]) - assert( abs(mdot(g_a,[g_aa]) - sdot(g_a,g_aa)).max() == 0) - Z = g_aa - V_aa = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('a', 'e') - - K_ae = + mdot(f_2,[V_a,V_e]) - L_ae = + mdot(g_aa,[g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ae) + K_ae - g_ae = - sdot(A_inv, const) - - if order < max_order: - Y = L_ae + mdot(g_a,[g_ae]) - Z = g_ae - V_ae = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('e', 'e') - - K_ee = + mdot(f_2,[V_e,V_e]) - L_ee = + mdot(g_aa,[g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ee) + K_ee - g_ee = - sdot(A_inv, const) - - if order < max_order: - Y = L_ee + mdot(g_a,[g_ee]) - Z = g_ee - V_ee = build_V(Y,Z,(n_a,n_e)) - - # manual - I = np.eye(n_v,n_v) - M_inv = np.linalg.inv( sdot(f1_A,g_a+I) + f1_B ) - K_ss = mdot(f_2[:,:n_v,:n_v],[g_e,g_e]) + sdot( f1_A, g_ee ) - rhs = - np.tensordot( K_ss, Sigma, axes=((1,2),(0,1)) ) #- mdot(h_2,[V_s,V_s]) - if derivatives_ss: - f_ss = derivatives_ss[0] - rhs -= f_ss - g_ss = sdot(M_inv,rhs) - ghs2 = g_ss/2 - - - - if max_order == 2: - d = { - 'ev': ev, - 'g_a': g_a, - 'g_e': g_e, - 'g_aa': g_aa, - 'g_ae': g_ae, - 'g_ee': g_ee, - 'g_ss': g_ss - } - return d - # /manual - - #----------Computing order 3 - - order = 3 - - #--- Computing derivatives ('a', 'a', 'a') - K_aaa = + 3*mdot(f_2,[V_a,V_aa]) + mdot(f_3,[V_a,V_a,V_a]) - L_aaa = + 3*mdot(g_aa,[g_a,g_aa]) - - #K_aaa = 2*( mdot(f_2,[V_aa,V_a]) ) + mdot(f_2,[V_a,V_aa]) + mdot(f_3,[V_a,V_a,V_a]) - #L_aaa = 2*( mdot(g_aa,[g_aa,g_a]) ) + mdot(g_aa,[g_a,g_aa]) - #K_aaa = ( mdot(f_2,[V_aa,V_a]) + mdot(f_2,[V_a,V_aa]) )*3.0/2.0 + mdot(f_3,[V_a,V_a,V_a]) - #L_aaa = ( mdot(g_aa,[g_aa,g_a]) + mdot(g_aa,[g_a,g_aa]) )*3.0/2.0 - - - #K_aaa = (K_aaa + K_aaa.swapaxes(3,2) + K_aaa.swapaxes(1,2) + K_aaa.swapaxes(1,2).swapaxes(2,3) + K_aaa.swapaxes(1,3) + K_aaa.swapaxes(1,3).swapaxes(2,3) )/6 - #L_aaa = (L_aaa + L_aaa.swapaxes(3,2) + L_aaa.swapaxes(1,2) + L_aaa.swapaxes(1,2).swapaxes(2,3) + L_aaa.swapaxes(1,3) + L_aaa.swapaxes(1,3).swapaxes(2,3) )/6 - - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aaa + sdot(f_d,L_aaa) - - - - if mlab == None: - g_aaa = solve_sylvester(A,B,C,D) - # this is much much faster - else: - n_d = D.ndim - 1 - n_v = C.shape[1] - CC = np.kron(np.kron(C,C),C) - DD = D.reshape( n_v, n_v**n_d ) - [err,E] = mlab.gensylv(3,A,B,C,DD,nout=2) - g_aaa = E.reshape((n_v,n_v,n_v,n_v)) - - #res = sdot(A,g_aaa) + sdot(B, mdot(g_aaa,[C,C,C])) - D - #print 'res : ' + str( abs(res).max() ) - - - if order < max_order: - Y = L_aaa + mdot(g_a,[g_aaa]) + mdot(g_aaa,[g_a,g_a,g_a]) - Z = g_aaa - V_aaa = build_V(Y,Z,(n_a,n_e)) - - # we transform g_aaa into a symmetric multilinear form - g_aaa = (g_aaa + g_aaa.swapaxes(3,2) + g_aaa.swapaxes(1,2) + g_aaa.swapaxes(1,2).swapaxes(2,3) + g_aaa.swapaxes(1,3) + g_aaa.swapaxes(1,3).swapaxes(2,3) )/6 - - #--- Computing derivatives ('a', 'a', 'e') - - K_aae = + mdot(f_2,[V_aa,V_e]) + 2*mdot(f_2,[V_a,V_ae]) + mdot(f_3,[V_a,V_a,V_e]) - L_aae = + mdot(g_aa,[g_aa,g_e]) + 2*mdot(g_aa,[g_a,g_ae]) + mdot(g_aaa,[g_a,g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aae) + K_aae - g_aae = - sdot(A_inv, const) - - if order < max_order: - Y = L_aae + mdot(g_a,[g_aae]) - Z = g_aae - V_aae = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('a', 'e', 'e') - - K_aee = + 2*mdot(f_2,[V_ae,V_e]) + mdot(f_2,[V_a,V_ee]) + mdot(f_3,[V_a,V_e,V_e]) - L_aee = + 2*mdot(g_aa,[g_ae,g_e]) + mdot(g_aa,[g_a,g_ee]) + mdot(g_aaa,[g_a,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aee) + K_aee - g_aee = - sdot(A_inv, const) - - if order < max_order: - Y = L_aee + mdot(g_a,[g_aee]) - Z = g_aee - V_aee = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('e', 'e', 'e') - - K_eee = + 3*mdot(f_2,[V_e,V_ee]) + mdot(f_3,[V_e,V_e,V_e]) - L_eee = + 3*mdot(g_aa,[g_e,g_ee]) + mdot(g_aaa,[g_e,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_eee) + K_eee - g_eee = - sdot(A_inv, const) - - if order < max_order: - Y = L_eee + mdot(g_a,[g_eee]) - Z = g_eee - V_eee = build_V(Y,Z,(n_a,n_e)) - - - #################################### - ## Compute sigma^2 correction term # - #################################### - - # ( a s s ) - - A = f_a + sdot(f_d,g_a) - I_e = np.eye(n_e) - - Y = g_e - Z = np.zeros((n_a,n_e)) - V_s = build_V(Y,Z,(n_a,n_e)) - - Y = mdot( g_ae, [g_a, I_e] ) - Z = np.zeros((n_a,n_a,n_e)) - V_as = build_V(Y,Z,(n_a,n_e)) - - Y = sdot(g_a,g_ss) + g_ss + np.tensordot(g_ee,Sigma) - Z = g_ss - V_ss = build_V(Y,Z,(n_a,n_e)) - - K_ass_1 = 2*mdot(f_2,[V_as,V_s] ) + mdot(f_3,[V_a,V_s,V_s]) - K_ass_1 = np.tensordot(K_ass_1,Sigma) - - K_ass_2 = mdot( f_2, [V_a,V_ss] ) - - K_ass = K_ass_1 + K_ass_2 - - L_ass = mdot( g_aa, [g_a, g_ss]) + np.tensordot( mdot(g_aee,[g_a, I_e, I_e]), Sigma) - - D = K_ass + sdot(f_d,L_ass) - - if derivatives_ss: - f_1ss = derivatives_ss[1] - D += mdot( f_1ss, [V_a]) + sdot( f1_A, mdot( g_aa, [ g_a , g_ss ]) ) - - g_ass = solve_sylvester( A, B, C, D) - - - # ( e s s ) - - A = f_a + sdot(f_d,g_a) - A_inv = np.linalg.inv(A) - I_e = np.eye(n_e) - - Y = g_e - Z = np.zeros((n_a,n_e)) - V_s = build_V(Y,Z,(n_a,n_e)) - - Y = mdot( g_ae, [g_e, I_e] ) - Z = np.zeros((n_a,n_e,n_e)) - V_es = build_V(Y,Z,(n_a,n_e)) - - Y = sdot(g_a,g_ss) + g_ss + np.tensordot(g_ee,Sigma) - Z = g_ss - V_ss = build_V(Y,Z,(n_a,n_e)) - - K_ess_1 = 2*mdot(f_2,[V_es,V_s] ) + mdot(f_3,[V_e,V_s,V_s]) - K_ess_1 = np.tensordot(K_ess_1,Sigma) - - K_ess_2 = mdot( f_2, [V_e,V_ss] ) - - K_ess = K_ess_1 + K_ess_2 - - L_ess = mdot( g_aa, [g_e, g_ss]) + np.tensordot( mdot(g_aee,[g_e, I_e, I_e]), Sigma) - L_ess += mdot( g_ass, [g_e]) - - D = K_ess + sdot(f_d,L_ess) - - g_ess = sdot( A_inv, -D) - - if max_order == 3: - d = {'ev':ev,'g_a':g_a,'g_e':g_e, 'g_aa':g_aa, 'g_ae':g_ae, 'g_ee':g_ee, - 'g_aaa':g_aaa, 'g_aae':g_aae, 'g_aee':g_aee, 'g_eee':g_eee, 'g_ss':g_ss, 'g_ass':g_ass,'g_ess':g_ess} - return d - - - -def new_solver_with_p(derivatives, sizes, max_order=2): - - if max_order == 1: - [f_0,f_1] = derivatives - elif max_order == 2: - [f_0,f_1,f_2] = derivatives - elif max_order == 3: - [f_0,f_1,f_2,f_3] = derivatives - derivs = derivatives - - f = derivs - #n = f[0].shape[0] # number of variables - #s = f[1].shape[1] - 3*n - [n_v,n_s,n_p] = sizes - - n = n_v - - f1_A = f[1][:,:n] - f1_B = f[1][:,n:(2*n)] - f1_C = f[1][:,(2*n):(3*n)] - f1_D = f[1][:,(3*n):((3*n)+n_s)] - f1_E = f[1][:,((3*n)+n_s):] - - ## first order - [ev,g_x] = second_order_solver(f1_A,f1_B,f1_C) - - mm = np.dot(f1_A, g_x) + f1_B - - g_u = - np.linalg.solve( mm , f1_D ) - g_p = - np.linalg.solve( mm , f1_E ) - - d = { - 'ev':ev, - 'g_a':g_x, - 'g_e':g_u, - 'g_p':g_p - } - - if max_order == 1: - return d - - # we need it for higher order - - V_x = np.concatenate( [np.dot(g_x,g_x),g_x,np.eye(n_v),np.zeros((n_s,n_v)), np.zeros((n_p,n_v))] ) - V_u = np.concatenate( [np.dot(g_x,g_u),g_u,np.zeros((n_v,n_s)),np.eye(n_s), np.zeros((n_p,n_s))] ) - V_p = np.concatenate( [np.dot(g_x,g_p),g_p,np.zeros((n_v,n_p)),np.zeros((n_s,n_p)), np.eye(n_p)] ) - V = [None, [V_x,V_u]] - - # Translation - n_a = n_v - n_e = n_s - n_p = g_p.shape[1] - - f_1 = f[1] - f_2 = f[2] - f_d = f1_A - f_a = f1_B - f_h = f1_C - f_u = f1_D - V_a = V_x - V_e = V_u - g_a = g_x - g_e = g_u - - # Once for all ! - A = f_a + sdot(f_d,g_a) - B = f_d - C = g_a - A_inv = np.linalg.inv(A) - - #----------Computing order 2 - - order = 2 - - #--- Computing derivatives ('a', 'a') - - K_aa = + mdot(f_2,[V_a,V_a]) - L_aa = np.zeros((n_v, n_a, n_a)) - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aa + sdot(f_d,L_aa) - g_aa = solve_sylvester(A,B,C,D) - - if order < max_order: - Y = L_aa + mdot(g_a,[g_aa]) + mdot(g_aa,[g_a,g_a]) - Z = g_aa - V_aa = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'e') - - K_ae = + mdot(f_2,[V_a,V_e]) - L_ae = + mdot(g_aa,[g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ae) + K_ae - g_ae = - sdot(A_inv, const) - - if order < max_order: - Y = L_ae + mdot(g_a,[g_ae]) - Z = g_ae - V_ae = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'p') - - K_ap = + mdot(f_2,[V_a,V_p]) - L_ap = + mdot(g_aa,[g_a,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ap) + K_ap - g_ap = - sdot(A_inv, const) - - if order < max_order: - Y = L_ap + mdot(g_a,[g_ap]) - Z = g_ap - V_ap = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'e') - - K_ee = + mdot(f_2,[V_e,V_e]) - L_ee = + mdot(g_aa,[g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ee) + K_ee - g_ee = - sdot(A_inv, const) - - if order < max_order: - Y = L_ee + mdot(g_a,[g_ee]) - Z = g_ee - V_ee = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'p') - - K_ep = + mdot(f_2,[V_e,V_p]) - L_ep = + mdot(g_aa,[g_e,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ep) + K_ep - g_ep = - sdot(A_inv, const) - - if order < max_order: - Y = L_ep + mdot(g_a,[g_ep]) - Z = g_ep - V_ep = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('p', 'p') - - K_pp = + mdot(f_2,[V_p,V_p]) - L_pp = + mdot(g_aa,[g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_pp) + K_pp - g_pp = - sdot(A_inv, const) - - if order < max_order: - Y = L_pp + mdot(g_a,[g_pp]) - Z = g_pp - V_pp = build_V(Y,Z,(n_a,n_e,n_p)) - - - d.update({ - 'g_aa':g_aa, - 'g_ae':g_ae, - 'g_ee':g_ee, - 'g_ap':g_ap, - 'g_ep':g_ep, - 'g_pp':g_pp - }) - if max_order == 2: - return d - - #----------Computing order 3 - - order = 3 - - #--- Computing derivatives ('a', 'a', 'a') - - K_aaa = + 3*mdot(f_2,[V_a,V_aa]) + mdot(f_3,[V_a,V_a,V_a]) - L_aaa = + 3*mdot(g_aa,[g_a,g_aa]) - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aaa + sdot(f_d,L_aaa) - g_aaa = solve_sylvester(A,B,C,D) - - if order < max_order: - Y = L_aaa + mdot(g_a,[g_aaa]) + mdot(g_aaa,[g_a,g_a,g_a]) - Z = g_aaa - V_aaa = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'a', 'e') - - K_aae = + mdot(f_2,[V_aa,V_e]) + 2*mdot(f_2,[V_a,V_ae]) + mdot(f_3,[V_a,V_a,V_e]) - L_aae = + mdot(g_aa,[g_aa,g_e]) + 2*mdot(g_aa,[g_a,g_ae]) + mdot(g_aaa,[g_a,g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aae) + K_aae - g_aae = - sdot(A_inv, const) - - if order < max_order: - Y = L_aae + mdot(g_a,[g_aae]) - Z = g_aae - V_aae = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'a', 'p') - - K_aap = + mdot(f_2,[V_aa,V_p]) + 2*mdot(f_2,[V_a,V_ap]) + mdot(f_3,[V_a,V_a,V_p]) - L_aap = + mdot(g_aa,[g_aa,g_p]) + 2*mdot(g_aa,[g_a,g_ap]) + mdot(g_aaa,[g_a,g_a,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aap) + K_aap - g_aap = - sdot(A_inv, const) - - if order < max_order: - Y = L_aap + mdot(g_a,[g_aap]) - Z = g_aap - V_aap = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'e', 'e') - - K_aee = + 2*mdot(f_2,[V_ae,V_e]) + mdot(f_2,[V_a,V_ee]) + mdot(f_3,[V_a,V_e,V_e]) - L_aee = + 2*mdot(g_aa,[g_ae,g_e]) + mdot(g_aa,[g_a,g_ee]) + mdot(g_aaa,[g_a,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aee) + K_aee - g_aee = - sdot(A_inv, const) - - if order < max_order: - Y = L_aee + mdot(g_a,[g_aee]) - Z = g_aee - V_aee = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'e', 'p') - ll = [ mdot(f_2,[V_ae,V_p]) , mdot(f_2,[V_ap,V_e]), mdot(f_2,[V_a,V_ep]) , mdot(f_3,[V_a,V_e,V_p]) ] - l = [ mdot(f_2,[V_ae,V_p]) , mdot(f_2,[V_ap,V_e]).swapaxes(2,3) , mdot(f_2,[V_a,V_ep]) , mdot(f_3,[V_a,V_e,V_p]) ] - - K_aep = + mdot(f_2,[V_ae,V_p]) + mdot(f_2,[V_ap,V_e]).swapaxes(2,3) + mdot(f_2,[V_a,V_ep]) + mdot(f_3,[V_a,V_e,V_p]) - L_aep = + mdot(g_aa,[g_ae,g_p]) + mdot(g_aa,[g_ap,g_e]).swapaxes(2,3) + mdot(g_aa,[g_a,g_ep]) + mdot(g_aaa,[g_a,g_e,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aep) + K_aep - g_aep = - sdot(A_inv, const) - - if order < max_order: - Y = L_aep + mdot(g_a,[g_aep]) - Z = g_aep - V_aep = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'p', 'p') - - K_app = + 2*mdot(f_2,[V_ap,V_p]) + mdot(f_2,[V_a,V_pp]) + mdot(f_3,[V_a,V_p,V_p]) - L_app = + 2*mdot(g_aa,[g_ap,g_p]) + mdot(g_aa,[g_a,g_pp]) + mdot(g_aaa,[g_a,g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_app) + K_app - g_app = - sdot(A_inv, const) - - if order < max_order: - Y = L_app + mdot(g_a,[g_app]) - Z = g_app - V_app = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'e', 'e') - - K_eee = + 3*mdot(f_2,[V_e,V_ee]) + mdot(f_3,[V_e,V_e,V_e]) - L_eee = + 3*mdot(g_aa,[g_e,g_ee]) + mdot(g_aaa,[g_e,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_eee) + K_eee - g_eee = - sdot(A_inv, const) - - if order < max_order: - Y = L_eee + mdot(g_a,[g_eee]) - Z = g_eee - V_eee = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'e', 'p') - - K_eep = + mdot(f_2,[V_ee,V_p]) + 2*mdot(f_2,[V_e,V_ep]) + mdot(f_3,[V_e,V_e,V_p]) - L_eep = + mdot(g_aa,[g_ee,g_p]) + 2*mdot(g_aa,[g_e,g_ep]) + mdot(g_aaa,[g_e,g_e,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_eep) + K_eep - g_eep = - sdot(A_inv, const) - - if order < max_order: - Y = L_eep + mdot(g_a,[g_eep]) - Z = g_eep - V_eep = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'p', 'p') - - K_epp = + 2*mdot(f_2,[V_ep,V_p]) + mdot(f_2,[V_e,V_pp]) + mdot(f_3,[V_e,V_p,V_p]) - L_epp = + 2*mdot(g_aa,[g_ep,g_p]) + mdot(g_aa,[g_e,g_pp]) + mdot(g_aaa,[g_e,g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_epp) + K_epp - g_epp = - sdot(A_inv, const) - - if order < max_order: - Y = L_epp + mdot(g_a,[g_epp]) - Z = g_epp - V_epp = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('p', 'p', 'p') - - K_ppp = + 3*mdot(f_2,[V_p,V_pp]) + mdot(f_3,[V_p,V_p,V_p]) - L_ppp = + 3*mdot(g_aa,[g_p,g_pp]) + mdot(g_aaa,[g_p,g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ppp) + K_ppp - g_ppp = - sdot(A_inv, const) - - if order < max_order: - Y = L_ppp + mdot(g_a,[g_ppp]) - Z = g_ppp - V_ppp = build_V(Y,Z,(n_a,n_e,n_p)) - - - d.update({ - 'g_aaa':g_aaa, - 'g_aae':g_aae, - 'g_aee':g_aee, - 'g_eee':g_eee, - 'g_aap':g_aap, - 'g_aep':g_aep, - 'g_eep':g_eep, - 'g_app':g_app, - 'g_epp':g_epp, - 'g_ppp':g_ppp - }) - - return d - - -def sigma2_correction(d,order,n_v,n_s,f1_A,f1_B, f_2, h_2, f_ss, Sigma): - - resp = dict() - - if order >= 2: - g_a = d['g_a'] - g_e = d['g_e'] - g_aa = d['g_aa'] - g_ae = d['g_ae'] - g_ee = d['g_ee'] - I = np.eye(n_v,n_v) - A_inv = np.linalg.inv( sdot(f1_A,g_a+I) + f1_B ) - V_s = np.zeros( (n_v*3+n_s+1,)) - V_s[-1] = 1 - K_ss = mdot(f_2[:,:n_v,:n_v],[g_e,g_e]) + sdot( f1_A, g_ee ) - rhs = np.tensordot( K_ss,Sigma) + mdot(h_2,[V_s,V_s]) - sigma2 = sdot(A_inv,-rhs) / 2 - - ds2 = sigma2 - d['g_ss'] - resp['sigma2'] = sigma2 - - if order == 3: - A = sdot(f1_A,g_a) + f1_B - A_inv = np.linalg.inv( A ) - B = f1_A - C = g_a - V_x = np.row_stack( [ - np.dot(g_a,g_a), - g_a, - I, - np.zeros((n_s,n_v)) - ]) - - D = mdot( f_ss, [V_x]) + sdot( f1_A, mdot( g_aa, [ g_a , ds2 ]) ) - dsigma2 = solve_sylvester(A,B,C,D) - resp['dsigma2'] = dsigma2 - resp['D'] = D - resp['f_ss'] = f_ss - return resp - -def build_V(X,Y,others): - - ddim = X.shape[1:] - new_dims = [(s,) + ddim for s in others] - l = [X,Y] + [np.zeros(d) for d in new_dims] - return np.concatenate(l, axis=0) - - - -if __name__ == '__main__': - - from dolo import * - - fname = "/home/pablo/Programming/econforge/dolo/examples/models/compat/rbc_dynare.yaml" - - - smodel = yaml_import(fname, return_symbolic=True) - - infos = { - 'name' : 'anonymous', - 'type' : 'dynare' - } - from dolo.compiler.model_dynare import NumericModel - - model = NumericModel(smodel, infos=infos) - - dr = solve_decision_rule(model, order=3) - - print(dr) diff --git a/trash/dolo/algos/dynare/simulations.py b/trash/dolo/algos/dynare/simulations.py deleted file mode 100644 index da483ed5..00000000 --- a/trash/dolo/algos/dynare/simulations.py +++ /dev/null @@ -1,143 +0,0 @@ -import numpy as np - -# from dolo.numeric.decision_rules import * - -# def impulse_response_function(decision_rule, shock, variables = None, horizon=40, order=1, output='deviations', plot=False): -# -# model = decision_rule.model -# -# -# if order > 1: -# raise Exception('irfs, for order > 1 not implemented') -# -# dr = decision_rule -# A = dr['g_a'] -# B = dr['g_e'] -# Sigma = dr['Sigma'] -# -# [n_v, n_s] = [ len(model.symbols['variables']), len(model.symbols['shocks']) ] -# -# -# if isinstance(shock,int): -# i_s = shock -# elif isinstance(shock,str): -# i_s = model.symbols['shocks'].index( shock ) -# else: -# raise Exception("Unknown type for shock") -# -# E0 = np.zeros( n_s ) -# E0[i_s] = np.sqrt(dr['Sigma'][i_s,i_s]) -# -# E = E0*0 -# -# # RSS = dr.risky_ss() -# RSS = dr['ys'] -# -# simul = np.zeros( (n_v, horizon+1) ) -# -# start = dr( RSS, E0 ) -# simul[:,0] = start -# for i in range(horizon): -# simul[:,i+1] = dr( simul[:,i], E ) -# -# # TODO: change the correction so that it corresponds to the risky steady-state -# constant = np.tile(RSS, ( horizon+1, 1) ).T -# if output == 'deviations': -# simul = simul - constant -# elif output == 'logs': -# simul = np.log((simul-constant)/constant) -# elif output == 'percentages': -# simul = (simul-constant)/constant*100 -# elif output == 'levels': -# pass -# else: -# raise Exception("Unknown output type") -# -# -# if variables: -# variables = [(str(v)) for v in variables] -# ind_vars = [model.symbols['variables'].index( v ) for v in variables] -# simul = simul[ind_vars, :] -# else: -# variables = model.symbols['variables'] -# -# -# x = np.linspace(0,horizon,horizon+1) -# -# if plot: -# from matplotlib import pylab -# pylab.clf() -# # pylab.title('Impulse-Response Function for ${var}$'.format(var=shock.__latex__())) -# pylab.title('Impulse-Response Function for ${var}$'.format(var=shock)) -# for k in range(len(variables)): -# # pylab.plot(x[1:], simul[k,1:],label='$'+variables[k]._latex_()+'$' ) -# pylab.plot(x[1:], simul[k,1:],label=variables[k] ) -# pylab.plot(x,x*0,'--',color='black') -# pylab.xlabel('$t$') -# if output == 'percentages': -# pylab.ylabel('% deviations from the steady-state') -# elif output == 'deviations': -# pylab.ylabel('Deviations from the steady-state') -# elif output == 'levels': -# pylab.ylabel('Levels') -# pylab.legend() -# # if dolo.config.save_plots: -# # filename = 'irf_' + str(shock) + '__' + '_' + str.join('_',[str(v) for v in variables]) -# # pylab.savefig(filename) # not good... -# # else: -# pylab.show() -# -# import pandas -# sim = pandas.DataFrame(simul.T, columns=variables, index=range(horizon+1)) -# return sim - -def simulate(decision_rule, horizon=40, start=None, shock=None, seed=None, n_exp=0): - - dr = decision_rule - model = dr.model - - if n_exp==0: - n_exp = 1 - irf = True - else: - irf = False - - [n_v, n_s] = [ len(model.symbols['variables']), len(model.symbols['shocks']) ] - - Sigma = dr['Sigma'] - if seed: - np.random.seed(seed) - - if irf: - E = np.zeros((n_exp,horizon+1,n_s)) - else: - E = np.random.multivariate_normal((0,)*n_s,Sigma,(n_exp,horizon+1)) - E[:,0,:] = 0 - - if shock is not None: - E[:,1,:] = np.array(shock).ravel()[None,:] - - simul = np.zeros( (n_exp,horizon+1, n_v) ) - - if start in ('risky', None): - RSS = dr.risky_ss() - start = RSS - elif start == 'deterministic': - start = dr['ys'] - - # TODO make multiple simulations more efficient - for n in range(n_exp): - simul[n,0,:] = start - for i in range(horizon): - simul[n,i+1,:] = dr( simul[n,i,:], E[n,i+1,:] ) - - all_vars = model.symbols['variables'] + model.symbols['shocks'] - sims = np.concatenate([simul, E], axis=2) - - import pandas - if n_exp == 1: - sim = pandas.DataFrame(sims[0,:,:], columns=all_vars, index=range(horizon+1)) - else: - sim = pandas.Panel(sims, items=range(n_exp), major_axis=range(horizon+1), minor_axis=all_vars) - - return sim diff --git a/trash/dolo/algos/dynare/steady_state.py b/trash/dolo/algos/dynare/steady_state.py deleted file mode 100644 index b1368b80..00000000 --- a/trash/dolo/algos/dynare/steady_state.py +++ /dev/null @@ -1,9 +0,0 @@ -def residuals(model, calib=None): - - if calib is None: - calib = model.calibration - - y,e,p = model.calibration['variables', 'shocks', 'parameters'] - res = model.functions['f_static'](y, p) - - return {'dynare': res} diff --git a/trash/dolo/misc/commands.py b/trash/dolo/misc/commands.py deleted file mode 100644 index 18e29999..00000000 --- a/trash/dolo/misc/commands.py +++ /dev/null @@ -1,23 +0,0 @@ - - - -def dynare(filename, order=1, plot=True): - - from trash.dolo.misc.modfile import dynare_import - from trash.dolo.numeric.perturbations import solve_decision_rule - from dolo.numeric.decision_rules import stoch_simul - - model = dynare_import(filename) - - dr = solve_decision_rule - - simul = stoch_simul(dr, plot=plot) - - return dict( - model = model, - dr = dr, - simul = simul - ) - - - diff --git a/trash/dolo/misc/matlab.py b/trash/dolo/misc/matlab.py deleted file mode 100644 index 53b3a8c5..00000000 --- a/trash/dolo/misc/matlab.py +++ /dev/null @@ -1,117 +0,0 @@ -import sympy - -from trash.dolo.symbolic.symbolic import Variable,Parameter - -def convert_struct_to_dict(s): - # imperfect but enough for now - if len(s.dtype) == 0: - if s.shape == (1,): - return str(s[0]) - elif s.shape == (0,0): - return [] - elif s.shape == (1,1): - return s[0,0] - else: - return s - else: - # we suppose that we found a structure - d = dict() - ss = s[0,0] # actual content of the structure - for n in ss.dtype.names: - d[n] = convert_struct_to_dict(ss[n]) - return d - -def send_to_matlab(model,interactive=True,rmtemp=False,append=""): - tempname = model.fname + ".mod" - f = file(tempname,'w') - modtext = model.export_to_modfile(options={},append=append) - f.write(modtext) - f.close() - main_file = "%% main file for model : %s\n\n" % (model.fname) - main_file += "dynare('%s')\n" % (model.fname) - main_file += "save -v6 %s M_ oo_ options_\n" % (model.fname + "_res.mat") - main_file += "exit" - - g = file( model.fname + "_main.m",'w') - g.write(main_file) - g.close() - - if interactive: - pass - #import os - #command = 'gnome-terminal -e "matlab --persist --nodesktop --nojvm --eval \\"%s\\""' % ( model.fname + "_main.m" ) - #os.system( command) - #return None - else: - import os - from scipy import io - command = "matlab -nodesktop -nojvm -r %s" % ( model.fname + "_main" ) - os.system(command) - res = io.loadmat( model.fname + "_res.mat") - return res - - -def send_to_dynare_2(solver,mlab,rmtemp=False): - tempname = solver.model.fname + ".mod" - f = file(tempname,'w') - modtext = solver.export_to_modfile(options={},append="") - f.write(modtext) - mlab.dynare(solver.model.fname) - f.close() - -def retrieve_results(mlab): - d_options = dict() - m_options = mlab._get('options_') - m_names = mlab.fieldnames(m_options) - for m_name in m_names: - #print(m_name) - name = str(mlab.cell2mat(m_name)) - print(name) - #name = m_name._[0] - #print(name) - #d_options[name] = getattr(m_options,name)#m_options.__getattr__(m_name) - return(d_options) - -def value_to_mat(v): - if isinstance(v,bool): - if v: - return '1' - else: - return '0' - elif ( isinstance(v,float) or isinstance(v,int) ): - return str(v) - elif isinstance(v,list): - l = [value_to_mat(i) for i in v] - classes = [i.__class__ for i in v] - if (str in classes) or (sympy.Symbol in classes) or (Parameter in classes) or (Variable in classes): #list contains at least one string - return 'strvcat(%s)' %str.join(' , ',l) - else: - return '[%s]' %str.join(' ; ',l) - elif isinstance(v,sympy.Matrix): - return '[%s]' %v.__repr__().replace('\n',';').replace(',',' ') - elif str(v.__class__) == "": - if len(v.shape) <= 2: - return str(v).replace('\n','').replace('] [',' ; ') - else: - return 'reshape( {0} , {1} )'.format( str(v.flatten('F')).replace('\n','') , str(v.shape).strip('()') ) - #raise Warning('list conversion to matlab not implemented (will be soon)') - else: - return "'%s'" %str(v) - - -def struct_to_mat(d,name = 'dict'): - # d must be a dictionary - txt = '%s = struct;\n' %name - for k in d.keys(): - if isinstance(k,float) or isinstance(k,int): - key = 'n%s' % str(k) - else: - key = str(k) - if isinstance(d[k],dict): - txt += struct_to_mat(d[k],key) - txt += "%s = setfield(%s,'%s',%s) ;\n" %(name, name , key, key) - #raise(Exception('recursive dictionaries not allowed yet')) - else: - v = value_to_mat(d[k]) - txt += "%s = setfield(%s,'%s',%s) ;\n" %(name, name , key, v) - return(txt) \ No newline at end of file diff --git a/trash/dolo/misc/modfile.py b/trash/dolo/misc/modfile.py deleted file mode 100644 index 7f158c2f..00000000 --- a/trash/dolo/misc/modfile.py +++ /dev/null @@ -1,320 +0,0 @@ -import re - -import sympy - -from dolo.symbolic.symbolic import Variable, Shock, Parameter, Equation - -line_regex = re.compile( - "(\s*)$|" # blank - +"\s*@#(.*)$|" # macro - +"\s*//[^#](.*)$|" # comment - +"\s*%[^#](.*)$|" # comment - +"\s*\[(.*)\]\s*$|" # tag - +"(\s*[^\s].*)$" # instruction -) -tag_regex = re.compile("\s*(\w+)\s*=\s*'(.*)'") - -def parse_dynare_text(txt,add_model=True,full_output=False, debug=False): - ''' - Imports the content of a modfile into the current interpreter scope - ''' - # here we call "instruction group", a string finishing by a semicolon - # an "instruction group" can have several lines - # a line can be - # - a comment //... - # - an old-style tag //$... - # - a new-style tag [key1='value1',..] - # - macro-instruction @#... - # A Modfile contains several blocks (in this order) : - # - an initblock defining variables, exovariables, parameters, initialization - # inside the initblock the order of declaration doesn't matter - # - a model block with two special lines (model; end;) - # - optional blocks (like endval, shocks) - # seperated by free matlab instructions in any order; - # - all other instructions are ignored - - otxt = txt - otxt = otxt.replace("\r\n","\n") - otxt = otxt.replace("^","**") - - # first, we remove end-of-line comments : they are definitely lost - regex = re.compile("(.+)//[^#](.*)") - def remove_end_comment(line): - res = regex.search(line) - if res: - l = res.groups(1)[0] - return(l) - else: - return line - txt = str.join("\n",map(remove_end_comment,otxt.split("\n"))) - - name_regex = re.compile("//\s*fname\s*=\s*'(.*)'") - m = name_regex.search(txt) - if m: - fname = m.group(1) - else: - fname = None - - - instruction_groups = [Instruction_group(s) for s in txt.split(";")] - - instructions = [ig.instruction for ig in instruction_groups] - - if debug: - print('Elementary instructions') - for i in instruction_groups: - print(i) - - try: - imodel = [re.compile('model(\(.*\)|)').match(e) is not None for e in instructions] - imodel = imodel.index(True) - #imodel = instructions.index("model") #this doesn't work for "MODEL" - iend = instructions.index("end") - model_block = instruction_groups[imodel:(iend+1)] - init_block = instruction_groups[0:imodel] - except: - raise Exception('Model block could not be found.') - - next_instructions = instructions[(iend+1):] - next_instruction_groups = instruction_groups[(iend+1):] - - if 'initval' in next_instructions: - iinitval = next_instructions.index('initval') - iend = next_instructions.index('end',iinitval) - matlab_block_1 = next_instruction_groups[0:iinitval] - initval_block = next_instruction_groups[iinitval:(iend+1)] - next_instruction_groups = next_instruction_groups[(iend+1):] - next_instructions = next_instructions[(iend+1):] - else: - initval_block = None - matlab_block_1 = None - - if 'endval' in next_instructions: - iendval = next_instructions.index('endval') - iend = next_instructions.index('end',iendval) - matlab_block_2 = next_instruction_groups[0:iendval] - endval_block = next_instruction_groups[iendval:(iend+1)] - next_instruction_groups = next_instruction_groups[(iend+1):] - next_instructions = next_instructions[(iend+1):] - else: - endval_block = None - matlab_block_2 = None - - # TODO : currently shocks block needs to follow initval, this restriction should be removed - if 'shocks' in next_instructions: - ishocks = next_instructions.index('shocks') - iend = next_instructions.index('end',ishocks) - matlab_block_3 = next_instruction_groups[0:ishocks] - shocks_block = next_instruction_groups[ishocks:(iend+1)] - next_instruction_groups = next_instruction_groups[(iend+1):] - next_instructions = next_instructions[(iend+1):] - else: - shocks_block = None - matlab_block_3 = None - - try: - init_regex = re.compile("(parameters |var |varexo |)(.*)") - var_names = [] - varexo_names = [] - parameters_names = [] - declarations = {} - for ig in init_block: - if ig.instruction != '': - m = init_regex.match(ig.instruction) - if not m: - raise Exception("Unexpected instruction in init block : " + str(ig.instruction)) - if m.group(1) == '': - [lhs,rhs] = m.group(2).split("=") - lhs = lhs.strip() - rhs = rhs.strip() - declarations[lhs] = rhs - else: - arg = m.group(2).replace(","," ") - names = [vn.strip() for vn in arg.split()] - if m.group(1).strip() == 'var': - dest = var_names - elif m.group(1).strip() == 'varexo': - dest = varexo_names - elif m.group(1).strip() == 'parameters': - dest = parameters_names - for n in names: - if not n in dest: - dest.append(n) - else: - raise Exception("symbol %s has already been defined".format(n)) - except Exception as e: - raise Exception('Init block could not be read : ' + str(e) ) - # the following instruction set the variables "variables","shocks","parameters" - - - variables = [] - for vn in var_names: - v = Variable(vn) - variables.append(v) - - shocks = [] - for vn in varexo_names: - s = Shock(vn) - shocks.append(s) - - parameters = [] - for vn in parameters_names: - p = Parameter(vn) - parameters.append(p) - - parse_dict = dict() - for v in variables + shocks + parameters: - parse_dict[v.name] = v - - special_symbols = [sympy.exp,sympy.log,sympy.sin,sympy.cos, sympy.atan, sympy.tan] - for s in special_symbols: - parse_dict[str(s)] = s - parse_dict['sqrt'] = sympy.sqrt - - - # Read parameters values - parameters_values = {} - for p in declarations: - try: - rhs = eval(declarations[p], parse_dict) - except Exception as e: - Exception("Impossible to evaluate parameter value : " + str(e)) - try: - lhs = eval(p,parse_dict) - except Exception as e: - # here we could declare p - raise e - parameters_values[lhs] = rhs - - - # Now we read the model block - model_tags = model_block[0].tags - equations = [] - for ig in model_block[1:-1]: - if ig.instruction != '': - teq = ig.instruction.replace('^',"**") - if '=' in teq: - teqlhs,teqrhs = teq.split("=") - else: - teqlhs = teq - teqrhs = '0' - eqlhs = eval(teqlhs, parse_dict) - eqrhs = eval(teqrhs, parse_dict) - eq = Equation(eqlhs,eqrhs) - eq.tags.update(ig.tags) - # if eq.tags.has_key('name'): - # eq.tags[] = ig.tags['name'] - equations.append(eq) - - # Now we read the initval block - init_values = {} - if initval_block != None: - for ig in initval_block[1:-1]: - if len(ig.instruction.strip()) >0: - try: - [lhs,rhs] = ig.instruction.split("=") - except Exception as e: - print(ig.instruction) - raise e - init_values[eval(lhs,parse_dict)] = eval(rhs,parse_dict) - - # Now we read the endval block - # I don't really care about the endval block ! - - end_values = {} - if endval_block != None: - for ig in endval_block[1:-1]: - [lhs,rhs] = ig.instruction.split("=") - end_values[eval(lhs)] = eval(rhs) - - # Now we read the shocks block - covariances = None - if shocks_block != None: - covariances = sympy.zeros(len(shocks)) - regex1 = re.compile("var (.*?),(.*?)=(.*)|var (.*?)=(.*)") - for ig in shocks_block[1:-1]: - m = regex1.match(ig.instruction) - if not m: - raise Exception("unrecognized instruction in block shocks : " + str(ig.instruction)) - if m.group(1) != None: - varname1 = m.group(1).strip() - varname2 = m.group(2).strip() - value = m.group(3).strip().replace("^","**") - elif m.group(4) != None: - varname1 = m.group(4).strip() - varname2 = varname1 - value = m.group(5).strip().replace("^","**") - i = varexo_names.index(varname1) - j = varexo_names.index(varname2) - covariances[i,j] = eval(value,parse_dict) - covariances[j,i] = eval(value,parse_dict) - - calibration = {} - calibration.update(parameters_values) - calibration.update(init_values) - symbols = {'variables': variables, 'shocks': shocks, 'parameters': parameters} - - from trash.dolo.symbolic.model import SModel - model = SModel({'dynare_block': equations}, symbols, calibration, covariances) - return model - - -class Instruction_group(): - - def __init__(self,s): - self.string = s - self.span = 1 + s.count("\n") - self.tags = {} - self.instruction = "" - self.instruction_type = None - self.process_lines() - return None - - def __repr__(self): - return str([self.instruction_type, self.instruction, self.tags]) - - def process_lines(self): - lines = [] - entire_instruction = "" - for l in self.string.split("\n"): - g = line_regex.match(l).groups() - matches = [i for i in range(len(g)) if g[i]!=None] - if len(matches) !=1 : - raise Exception( "Parsing error" ) - else: - i = matches[0] - self.instruction_type = i - if i == 0: - # this line is blank, just ignore it - lines.append(["blank",g[i]]) - if i == 1: - # this line contains macro instructions - lines.append(["macro",g[i]]) - if i in (2,3): - # there are regular comments on this line - lines.append(["comment",g[i]]) - if i == 4: - # new style tags - lines.append(["tags",g[i]]) - pairs = [tag_regex.match(kt).groups() for kt in g[i].split(',')] - for p in pairs: - if len(p)>1: - self.tags[p[0]] = p[1].strip() - if i == 5: - entire_instruction += (" " + g[i]) - lines.append(("instruction",g[i])) - self.instruction = entire_instruction.strip() - self.lines = lines - - - -def dynare_import(filename,full_output=False, debug=False): - '''Imports model defined in specified file''' - import os - basename = os.path.basename(filename) - fname = re.compile('(.*)\.(.*)').match(basename).group(1) - f = open(filename) - txt = f.read() - model = parse_dynare_text(txt,full_output=full_output, debug=debug) - model.name = fname - return model diff --git a/trash/dolo/misc/printing.py b/trash/dolo/misc/printing.py deleted file mode 100644 index 34fe9b76..00000000 --- a/trash/dolo/misc/printing.py +++ /dev/null @@ -1,104 +0,0 @@ -from sympy.printing.str import StrPrinter - - -def sympy_to_dynare_string(sexpr): - s = str(sexpr) - s = s.replace("==","=") - s = s.replace("**","^") - return(s) - - -class DoloPrinter(StrPrinter): - def _print_TSymbol(self,expr): - return expr.__str__() - -dp = DoloPrinter() - - - -############################################################################# -# The following functions compute the HTML representation of common objects # -############################################################################# - - - -def print_table( tab, col_names=None, row_names=None, header=False): - #table_style = "background-color: #A8FFBE;" - #td_style = "background-color: #F8FFBD; border: medium solid white; padding: 5px; text-align: right" - txt_lines = '' - for l in tab: - txt_columns = ['\t\t{0} '.format(str(c)) for c in l] - txt_columns = str.join('\n',txt_columns) - txt_lines += '\t\n{0}\n\n'.format( txt_columns ) - txt = '{0}
'.format(txt_lines) - return txt - -def print_array( obj,row_names=None,col_names=None): - import numpy - tab = numpy.atleast_2d( obj ) - resp = [[ "%.4f" %tab[i,j] for j in range(tab.shape[1]) ] for i in range(tab.shape[0]) ] - if row_names: - resp = [ [row_names[i]] + resp[i] for i in range(tab.shape[0]) ] - if col_names: - if row_names: - resp = [[''] +col_names] + resp - else: - resp = [col_names] + resp - return print_table(resp) - -def print_model( model, print_residuals=True): - - from sympy import latex - if print_residuals: - from trash.dolo.symbolic.model import compute_residuals - res = compute_residuals(model) - if len( model.equations_groups ) > 0: - if print_residuals: - eqs = [ ['', 'Equations','Residuals'] ] - else: - eqs = [ ['', 'Equations'] ] - for groupname in model.equations_groups: - eqg = model.equations_groups - eqs.append( [ groupname ,''] ) - if print_residuals: - eqs.extend([ ['','${}$'.format(latex(eq)),str(res[groupname][i])] for i,eq in enumerate(eqg[groupname]) ]) - else: - eqs.extend([ ['','${}$'.format(latex(eq))] for eq in eqg[groupname] ]) - txt = print_table( eqs, header = True) - return txt - - else: - if print_residuals: - table = [ (i+1, '${}$'.format(latex(eq)), str(res[i]) ) for i,eq in enumerate(model.equations)] - else: - table = [ (i+1, '${}$'.format(latex(eq)) ) for i,eq in enumerate(model.equations)] - txt = print_table([['','Equations']] + table, header=True) - return txt - -def print_cmodel( cmodel, print_residuals=True): - - if cmodel.model is None: - return str(cmodel) - - model = cmodel.model - - txt = 'Compiled model : {}\n'.format(model.name) - txt += print_model(model, print_residuals=print_residuals) - return txt - -#from pandas import - -def print_dynare_decision_rule( dr ): - - col_names = [v(-1) for v in dr.model.variables] + dr.model.shocks - row_names = dr.model.variables - - #first_order_coeffs = df(column_stack([dr['g_a'],dr['g_e']]), columns = col_names, index = row_names) - from numpy import column_stack - [inds, state_vars] = zip( *[ [i,v] for i,v in enumerate(dr.model.variables) if v in dr.model.state_variables] ) - first_order_coeffs = column_stack([dr['g_a'][:,inds],dr['g_e']]) - - col_names = [v(-1) for v in state_vars] + dr.model.shocks - row_names = dr.model.variables - txt = print_array(first_order_coeffs, col_names = col_names, row_names = row_names) - return txt \ No newline at end of file diff --git a/trash/dolo/misc/symbolic_interactive.py b/trash/dolo/misc/symbolic_interactive.py deleted file mode 100644 index 507000aa..00000000 --- a/trash/dolo/misc/symbolic_interactive.py +++ /dev/null @@ -1,343 +0,0 @@ -import inspect -import re - -from trash.dolo.symbolic.symbolic import Variable,Parameter,Shock,IndexedSymbol - -#def set_variables(s,names_dict={}): -# """ -# Creates symbolic variable with the name *s*. -# -- a string, either a single variable name, or -# a space separated list of variable names, or -# a list of variable names. -# NOTE: The new variable is both returned and automatically injected into -# the parent's *global* namespace. It's recommended not to use "var" in -# library code, it is better to use symbols() instead. -# EXEMPLES: -# """ -# -# frame = inspect.currentframe().f_back -# try: -# if not isinstance(s, list): -# s = re.split('\s|,', s) -# res = [] -# for t in s: -# # skip empty stringG -# if not t: -# continue -# if t in names_dict: -# latex_name = names_dict[t] -# else: -# latex_name = None -# sym = Variable(t,0,latex_name) -# frame.f_globals[t] = sym -# res.append(sym) -# res = list(res) -# if len(res) == 0: # var('') -# res = [] -# # otherwise var('a b ...') -# frame.f_globals['variables'] = res -# return res -# finally: -# del frame -# -#def add_variables(s,latex_names=None): -# """ -# The same as set_variables but doesn't replace the existing variables. -# """ -# frame = inspect.currentframe().f_back -# try: -# if not isinstance(s, list): -# s = re.split('\s|,', s) -# if latex_names <> None: -# sl = re.split(' ', latex_names) -# if len(sl)<> len(s): -# raise Exception, "You should supply one latex name per variable" -# res = [] -# for i in range(len(s)): -# t=s[i] -# # skip empty stringG -# if not t: -# continue -# if latex_names == None: -# sym = Variable(t,0) -# else: -# sym = Variable(t,0,latex_name=sl[i]) -# frame.f_globals[t] = sym -# res.append(sym) -# res = list(res) -# if len(res) == 0: # var('') -# res = [] -# # otherwise var('a b ...') -# frame.f_globals['variables'] += res -# return res -# finally: -# del frame -# -#def set_shocks(s,latex_names=None,names_dict={}): -# """ -# Creates symbolic variable with the name *s*. -# -- a string, either a single variable name, or -# a space separated list of variable names, or -# a list of variable names. -# NOTE: The new variable is both returned and automatically injected into -# the parent's *global* namespace. It's recommended not to use "var" in -# library code, it is better to use symbols() instead. -# EXAMPLES: -# """ -# -# frame = inspect.currentframe().f_back -# try: -# if not isinstance(s, list): -# s = re.split('\s|,', s) -# if latex_names <> None: -# sl = re.split(' ', latex_names) -# if len(sl)<> len(s): -# raise Exception, "You should supply one latex name per variable" -# res = [] -# for i in range(len(s)): -# t = s[i] -# # skip empty stringG -# if not t: -# continue -# if latex_names != None: -# sym = Shock(t,0,latex_name=sl[i]) -# elif t in names_dict: -# sym = Shock(t,0,latex_name=names_dict[t]) -# else: -# sym = Shock(t,0) -# frame.f_globals[t] = sym -# res.append(sym) -# res = list(res) -# if len(res) == 0: # var('') -# res = [] -# # otherwise var('a b ...') -# frame.f_globals['shocks'] = res -# return res -# finally: -# del frame -# -#def add_shocks(s,latex_names=None): -# """ -# The same as set_shocks but doesn't replace the existing variables. -# """ -# frame = inspect.currentframe().f_back -# try: -# if not isinstance(s, list): -# s = re.split('\s|,', s) -# if latex_names <> None: -# sl = re.split(' ', latex_names) -# if len(sl)<> len(s): -# raise Exception, "You should supply one latex name per variable" -# res = [] -# for i in range(len(s)): -# t=s[i] -# # skip empty stringG -# if not t: -# continue -# if latex_names == None: -# sym = Shock(t,0) -# else: -# sym = Shock(t,0,latex_name=sl[i]) -# frame.f_globals[t] = sym -# res.append(sym) -# res = list(res) -# if len(res) == 0: # var('') -# res = [] -# # otherwise var('a b ...') -# frame.f_globals['shocks'] += res -# return res -# finally: -# del frame -# -#def set_parameters(s,names_dict={}): -# """Create S symbolic variable with the name *s*. -# -- a string, either a single variable name, or -# a space separated list of variable names, or -# a list of variable names. -# NOTE: The new variable is both returned and automatically injected into -# the parent's *global* namespace. It's recommended not to use "var" in -# library code, it is better to use symbols() instead. -# EXAMPLES: -# """ -# -# frame = inspect.currentframe().f_back -# try: -# if not isinstance(s, list): -# s = re.split('\s|,', s) -# res = [] -# for t in s: -# # skip empty stringG -# if not t: -# continue -# if t in names_dict: -# latex_name = names_dict[t] -# else: -# latex_name = None -# sym = Parameter(t,latex_name) -# frame.f_globals[t] = sym -# res.append(sym) -# res = list(res) -# if len(res) == 0: # var('') -# res = [] -# # otherwise var('a b ...') -# frame.f_globals['parameters'] = res -# return res -# finally: -# del frame -# -#def add_parameters(s,latex_names=None): -# """ -# The same as set_variables but doesn't replace the existing variables. -# """ -# -# frame = inspect.currentframe().f_back -# try: -# if not isinstance(s, list): -# s = re.split('\s|,', s) -# if latex_names <> None: -# sl = re.split('\s|,', latex_names) -# if len(sl)<> len(s): -# raise Exception, "You should supply one latex name per variable" -# res = [] -# for i in range(len(s)): -# t=s[i] -# # skip empty stringG -# if not t: -# continue -# if latex_names == None: -# sym = Parameter(t) -# else: -# sym = Parameter(t,latex_name=sl[i]) -# frame.f_globals[t] = sym -# res.append(sym) -# res = list(res) -# if len(res) == 0: # var('') -# res = [] -# # otherwise var('a b ...') -# if frame.f_globals.get('parameters'): -# frame.f_globals['parameters'].extend(res) -# else: -# frame.f_globals['parameters'] = res -# return res -# finally: -# del frame - -#### new style ####### - -def def_variables(s): - """ - blabla - """ - - frame = inspect.currentframe().f_back - try: - if isinstance(s,str): - s = re.split('\s|,', s) - res = [] - for t in s: - # skip empty stringG - if not t: - continue - if t.count("@") > 0: - sym = IndexedSymbol(t,Variable) - t = t.strip('@') - else: - sym = Variable(t) - frame.f_globals[t] = sym - res.append(sym) - if frame.f_globals.get('variables_order'): - # we should avoid to declare symbols twice ! - frame.f_globals['variables_order'].extend(res) - else: - frame.f_globals['variables_order'] = res - return res - finally: - del frame - -def def_shocks(s): - """ - blabla - """ - - frame = inspect.currentframe().f_back - try: - if isinstance(s,str): - s = re.split('\s|,', s) - res = [] - for t in s: - # skip empty stringG - if not t: - continue - if t.count("@") > 0: - sym = IndexedSymbol(t,Shock) - t = t.strip('@') - else: - sym = Shock(t) - frame.f_globals[t] = sym - res.append(sym) - if frame.f_globals.get('shocks_order'): - # we should avoid to declare symbols twice ! - frame.f_globals['shocks_order'].extend(res) - else: - frame.f_globals['shocks_order'] = res - return res - finally: - del frame - - -def def_parameters(s): - """ - blabla - """ - - frame = inspect.currentframe().f_back - try: - if isinstance(s,str): - s = re.split('\s|,', s) - res = [] - for t in s: - # skip empty stringG - if not t: - continue - if t.count("@") > 0: - sym = IndexedSymbol(t,Parameter) - t = t.strip('@') - else: - sym = Parameter(t) - frame.f_globals[t] = sym - res.append(sym) - if frame.f_globals.get('parameters_order'): - # we should avoid to declare symbols twice ! - frame.f_globals['parameters_order'].extend(res) - else: - frame.f_globals['parameters_order'] = res - return res - finally: - del frame - - -def clear_all(): - """ - Clears all parameters, variables, and shocks defined previously - """ - - frame = inspect.currentframe().f_back - try: - if frame.f_globals.get('variables_order'): - # we should avoid to declare symbols twice ! - del frame.f_globals['variables_order'] - if frame.f_globals.get('parameters_order'): - # we should avoid to declare symbols twice ! - del frame.f_globals['parameters_order'] - finally: - del frame - - -def inject_symbols(symbs): - frame = inspect.currentframe().f_back - try: - for s in symbs: - sn = s.name - frame.f_globals[sn] = s - finally: - del frame diff --git a/trash/dolo/numeric/decision_rules.py b/trash/dolo/numeric/decision_rules.py deleted file mode 100644 index 0a65a0fd..00000000 --- a/trash/dolo/numeric/decision_rules.py +++ /dev/null @@ -1,504 +0,0 @@ -""" -This module contains classes representing decision rules -""" - -import numpy as np - -import dolo.config -from dolo.misc.caching import memoized - -class TaylorExpansion(dict): - - @property - @memoized - def order(self): - order = 0 - if self.get('g_a') is not None: - order = 1 - if self.get('g_aa') is not None: - order = 2 - if self.get('g_aaa') is not None: - order = 3 - return order - - - -class DynareDecisionRule(TaylorExpansion): - - - def __init__(self,d,model): - super(DynareDecisionRule,self).__init__(d) - self.model = model - self.dr_var_order_i = [model.variables.index(v) for v in self.dr_var_order] - self.dr_states_order = [v for v in self.dr_var_order if v in self.state_variables] - self.dr_states_i = [model.variables.index(v) for v in self.dr_states_order] - - @property - @memoized - def ghx(self): - ghx = self.get('g_a') - ghx = ghx[self.dr_var_order_i,:] - ghx = ghx[:,self.dr_states_i] - return ghx - - @property - @memoized - def ghu(self): - ghu = self.get('g_e') - ghu = ghu[self.dr_var_order_i,:] - return ghu - - @property - @memoized - def ghxx(self): - ghxx = self.get('g_aa') - ghxx = ghxx[self.dr_var_order_i,:] - ghxx = ghxx[:,self.dr_states_i,:] - ghxx = ghxx[:,:,self.dr_states_i] - n_v = ghxx.shape[0] - n_s = ghxx.shape[1] - return ghxx.reshape( (n_v,n_s*n_s) ) - - @property - @memoized - def ghxu(self): - ghxu = self.get('g_ae') - ghxu = ghxu[self.dr_var_order_i,:,:] - ghxu = ghxu[:,self.dr_states_i,:] - n_v = ghxu.shape[0] - n_s = ghxu.shape[1] - n_e = ghxu.shape[2] - return ghxu.reshape( (n_v,n_s*n_e )) - - @property - @memoized - def ghuu(self): - ghuu = self.get('g_ee') - ghuu = ghuu[self.dr_var_order_i,:] - n_v = ghuu.shape[0] - n_e = ghuu.shape[1] - return ghuu.reshape( (n_v,n_e*n_e) ) - - @property - @memoized - def ghs2(self): - ghs2 = self.get('g_ss') - ghs2 = ghs2[self.dr_var_order_i] - return ghs2 - - - @property - @memoized - def g_1(self): - - g_1x = self.ghx - g_1u = self.ghu - - if 'g_ass' in self: - correc_x_ss = self['g_ass'][self.dr_var_order_i,:][:,self.dr_states_i]/2 - correc_u_ss = self['g_ess'][self.dr_var_order_i,:]/2 - g_1x += correc_x_ss - g_1u += correc_u_ss - - return np.column_stack([g_1x,g_1u]) - - @property - @memoized - def g_2(self): - n_v = self['g_a'].shape[0] - n_s = len(self.dr_states_order) - n_e = self['g_e'].shape[1] - g_2 = np.zeros( (n_v,n_s+n_e,n_s+n_e) ) - - - g_2[:,n_s:,n_s:] = self.ghuu.reshape((n_v,n_e,n_e)) - g_2[:,:n_s,n_s:] = self.ghxu.reshape((n_v,n_s,n_e)) - g_2[:,n_s:,:n_s] = self.ghxu.reshape((n_v,n_s,n_e)).swapaxes(1,2) - g_2[:,:n_s,:n_s] = self.ghxx.reshape((n_v,n_s,n_s)) - return fold( g_2 ) / 2 - - @property - @memoized - def ghxxx(self): - ghxxx = self['g_aaa'][self.dr_var_order_i,...] - ghxxx = ghxxx[:,self.dr_states_i,:,:] - ghxxx = ghxxx[:,:,self.dr_states_i,:] - ghxxx = ghxxx[:,:,:,self.dr_states_i] - return ghxxx - - @property - @memoized - def ghxxu(self): - ghxxu = self['g_aae'][self.dr_var_order_i,...] - ghxxu = ghxxu[:,self.dr_states_i,:,:] - ghxxu = ghxxu[:,:,self.dr_states_i,:] - return ghxxu - - @property - @memoized - def ghxuu(self): - ghxuu = self['g_aee'][self.dr_var_order_i,...] - ghxuu = ghxuu[:,self.dr_states_i,:,:] - return ghxuu - - @property - @memoized - def ghuuu(self): - ghuuu = self['g_eee'][self.dr_var_order_i,...] - return ghuuu - - @property - @memoized - def g_3(self): - n_v = self['g_a'].shape[0] - n_s = len(self.dr_states_order) - n_e = self['g_e'].shape[1] - - ghxxx = self.ghxxx - ghxxu = self.ghxxu - ghxuu = self.ghxuu - ghuuu = self.ghuuu - - g_3 = np.zeros( (n_v,n_s+n_e,n_s+n_e,n_s+n_e)) - - g_3[:,:n_s,:n_s,:n_s] = ghxxx - - g_3[:,:n_s,:n_s,n_s:] = ghxxu - g_3[:,:n_s,n_s:,:n_s] = ghxxu.swapaxes(2,3) - g_3[:,n_s:,:n_s,:n_s] = ghxxu.swapaxes(1,3) - - g_3[:, :n_s, n_s:, n_s:] = ghxuu - g_3[:, n_s:, :n_s, n_s:] = ghxuu.swapaxes(1,2) - g_3[:, n_s:, n_s:, :n_s] = ghxuu.swapaxes(1,3) - - g_3[:,n_s:,n_s:,n_s:] = ghuuu - - g_3 = symmetrize(g_3) - - return fold( g_3 ) / 2 / 3 - - - def __str__(self): - txt = ''' -Decision rule (order {order}) : -{msg} - - States : {states} - - - Endogenous variables : {endo} - - - First order coefficients : - -{foc} -''' - mat = np.concatenate([self.ghx,self.ghu],axis=1) - if self.order > 1: - msg = '\n (Only first order derivatives are printed)\n' - else: - msg = '' - col_names = [str(v(-1)) for v in self.model.dr_var_order if v in self.model.state_variables] + [str(s) for s in self.model.shocks] - row_names = [str(v) for v in self.model.dr_var_order] - txt = txt.format( - msg = msg, - order=self.order, - states=str.join(' ', col_names), - endo=str.join(' ', row_names), - foc=mat - ) - return txt - - @memoized - def risky_ss(self): - oo = self.order - if oo == 1: - return self['ys'] - elif oo in (2,3): - #TODO: RSS for order 3 should be computed differently - import numpy.linalg - A = self['g_a'] - I = np.eye(A.shape[0]) - D = self['g_ss']/2 - dx = numpy.linalg.solve( I - A, D) - return self['ys'] + dx - - def gap_to_risky_steady_state(self,x): - from dolo.numeric.tensor import mdot - d = x - self['ys'] - res = self['ys'] - x - res += np.dot( self['g_a'],d ) - res += mdot(self['g_aa'],[d,d]) - res += mdot(self['g_aaa'],[d,d,d]) - res += self['g_ss']/2 - res += np.dot(self['g_ass'],d)/2 - return res - - def __call__(self,x,e,order=None): - from dolo.numeric.tensor import mdot - if order is None: - order = self.order - d = x - self['ys'] - res = self['ys'] + np.dot( self['g_a'], d ) - res += np.dot( self['g_e'], e ) - if ('g_aa' in self) or (order >=2): - res += mdot(self['g_aa'],[d,d])/2 - res += mdot(self['g_ae'],[d,e]) - res += mdot(self['g_ee'],[e,e])/2 - res += self['g_ss']/2 - if ('g_aaa' in self) or (order >=3): - res += mdot(self['g_aaa'],[d,d,d])/6 - res += mdot(self['g_aae'],[d,d,e])/3 - res += mdot(self['g_aee'],[d,e,e])/3 - res += mdot(self['g_eee'],[e,e,e])/6 - res += np.dot(self['g_ass'],d)/2 - res += np.dot(self['g_ess'],e)/2 - return res - - @property - def dyn_var_order(this): - self = this.model - # returns a list of dynamic variables ordered as in Dynare's dynamic function - if hasattr(self,'__dyn_var_order__') : - return self.__dyn_var_order__ - d = dict() - for eq in self.equations: - all_vars = eq.variables - for v in all_vars: - if not v.lag in d: - d[v.lag] = set() - d[v.lag].add(v) - maximum = max(d.keys()) - minimum = min(d.keys()) - ord = [] - for i in range(minimum,maximum+1): - if i in d.keys(): - ord += [v(i) for v in self.variables if v(i) in d[i]] - self.__dyn_var_order__ = ord - return ord - - @property - def state_variables(self): - return [v for v in self.model.variables if v(-1) in self.dyn_var_order ] - - @property - def dr_var_order(self): - dvo = self.dyn_var_order - purely_backward_vars = [v for v in self.model.variables if (v(1) not in dvo) and (v(-1) in dvo)] - purely_forward_vars = [v for v in self.model.variables if (v(-1) not in dvo) and (v(1) in dvo)] - static_vars = [v for v in self.model.variables if (v(-1) not in dvo) and (v(1) not in dvo) ] - mixed_vars = [v for v in self.model.variables if not v in purely_backward_vars+purely_forward_vars+static_vars ] - dr_order = static_vars + purely_backward_vars + mixed_vars + purely_forward_vars - return dr_order - -DecisionRule = DynareDecisionRule - -def theoretical_moments(dr,with_correlations=True): - maxit = 1000 - tol = 0.00001 - A = dr['g_a'] - B = dr['g_e'] - Sigma = dr['Sigma'] - M0 = np.dot(B,np.dot(Sigma,B.T)) - M1 = M0 - for i in range( maxit ): - M = M0 + np.dot( A, np.dot( M1, A.T ) ) - if abs( M - M1).max() < tol: - break - M1 = M - if not with_correlations: - return M - else: - cov = M - d = np.diag( 1/np.sqrt( np.diag(cov) ) ) - correl = np.dot(d, np.dot(cov,d.T) ) - return [M,correl] - - - - - -def symmetrize(tens): - return (tens + tens.swapaxes(3,2) + tens.swapaxes(1,2) + tens.swapaxes(1,2).swapaxes(2,3) + tens.swapaxes(1,3) + tens.swapaxes(1,3).swapaxes(2,3) )/6 - -def fold(tens): - from itertools import product - - if tens.ndim == 3: - n = tens.shape[0] - q = tens.shape[1] - assert( tens.shape[2] == q ) - - non_decreasing_pairs = [ (i,j) for (i,j) in product(range(q),range(q)) if i<=j ] - result = np.zeros( (n,len(non_decreasing_pairs) ) ) - for k,(i,j) in enumerate(non_decreasing_pairs): - result[:,k] = tens[:,i,j] - return result - - elif tens.ndim == 4: - n = tens.shape[0] - q = tens.shape[1] - assert( tens.shape[2] == q ) - assert( tens.shape[3] == q) - - non_decreasing_tuples = [ (i,j,k) for (i,j,k) in product(range(q),range(q),range(q)) if i<=j<=k ] - result = np.zeros( (n,len(non_decreasing_tuples) ) ) - for l,(i,j,k) in enumerate(non_decreasing_tuples): - result[:,l] = tens[:,i,j,k] - return result - - - -def impulse_response_function(decision_rule, shock, variables = None, horizon=40, order=1, output='deviations', plot=True): - - if order > 1: - raise Exception('irfs, for order > 1 not implemented') - - dr = decision_rule - A = dr['g_a'] - B = dr['g_e'] - Sigma = dr['Sigma'] - - [n_v, n_s] = [ len(dr.model.variables), len(dr.model.shocks) ] - - - if isinstance(shock,int): - i_s = shock - elif isinstance(shock,str): - from trash.dolo.symbolic.symbolic import Shock - shock = Shock(shock) - i_s = dr.model.shocks.index( shock ) - else: - i_s = shock - - E0 = np.zeros( n_s ) - E0[i_s] = np.sqrt(dr['Sigma'][i_s,i_s]) - - E = E0*0 - - RSS = dr.risky_ss() - - simul = np.zeros( (n_v, horizon+1) ) - - start = dr( RSS, E0 ) - simul[:,0] = start - for i in range(horizon): - simul[:,i+1] = dr( simul[:,i], E ) - - # TODO: change the correction so that it corresponds to the risky steady-state - constant = np.tile(RSS, ( horizon+1, 1) ).T - if output == 'deviations': - simul = simul - constant - elif output == 'logs': - simul = np.log((simul-constant)/constant) - elif output == 'percentages': - simul = (simul-constant)/constant*100 - elif output == 'levels': - pass - else: - raise Exception("Unknown output type") - - - if variables: - from trash.dolo.symbolic.symbolic import Variable - if not isinstance(variables,list): - variables = [variables] - variables = [Variable(str(v)) for v in variables] - ind_vars = [dr.model.variables.index( v ) for v in variables] - simul = simul[ind_vars, :] - - x = np.linspace(0,horizon,horizon+1) - - if plot: - from matplotlib import pylab - pylab.clf() - pylab.title('Impulse-Response Function for ${var}$'.format(var=shock.__latex__())) - for k in range(len(variables)): - pylab.plot(x[1:], simul[k,1:],label='$'+variables[k]._latex_()+'$' ) - pylab.plot(x,x*0,'--',color='black') - pylab.xlabel('$t$') - if output == 'percentages': - pylab.ylabel('% deviations from the steady-state') - elif output == 'deviations': - pylab.ylabel('Deviations from the steady-state') - elif output == 'levels': - pylab.ylabel('Levels') - pylab.legend() - if dolo.config.save_plots: - filename = 'irf_' + str(shock) + '__' + '_' + str.join('_',[str(v) for v in variables]) - pylab.savefig(filename) # not good... - else: - pylab.show() - - return simul - -def stoch_simul(decision_rule, variables = None, horizon=40, order=None, start=None, output='deviations', plot=True, seed=None): - -# if order > 1: -# raise Exception('irfs, for order > 1 not implemented') - - dr = decision_rule - - [n_v, n_s] = [ len(dr.model.variables), len(dr.model.shocks) ] - - Sigma = dr['Sigma'] - if seed: - np.random.seed(seed) - - E = np.random.multivariate_normal((0,)*n_s,Sigma,horizon) - E = E.T - - simul = np.zeros( (n_v, horizon+1) ) - RSS = dr.risky_ss() - if start is None: - start = RSS - - if not order: - order = dr.order - - simul[:,0] = start - for i in range(horizon): - simul[:,i+1] = dr( simul[:,i], E[:,i] ) - - # TODO: change the correction so that it corresponds to the risky steady-state - constant = np.tile(RSS, ( horizon+1, 1) ).T - if output == 'deviations': - simul = simul - constant - elif output == 'logs': - simul = np.log((simul-constant)/constant) - elif output == 'percentages': - simul = (simul-constant)/constant*100 - elif output == 'levels': - pass - else: - raise Exception("Unknown output type") - - if variables: - from trash.dolo.symbolic.symbolic import Variable - if not isinstance(variables,list): - variables = [variables] - variables = [Variable(str(v)) for v in variables] - ind_vars = [dr.model.variables.index( v ) for v in variables] - simul = simul[ind_vars, :] - - x = np.linspace(0,horizon,horizon+1) - - if plot: - from matplotlib import pylab - pylab.clf() - pylab.title('Stochastic simulation') - for k in range(len(variables)): - pylab.plot(x[1:], simul[k,1:],label='$'+variables[k]._latex_()+'$' ) - pylab.plot(x,x*0,'--',color='black') - pylab.xlabel('$t$') - if output == 'percentages': - pylab.ylabel('% deviations from the steady-state') - elif output == 'deviations': - pylab.ylabel('Deviations from the steady-state') - elif output == 'levels': - pylab.ylabel('Levels') - pylab.legend() - if dolo.config.save_plots: - filename = 'simul_' + '_' + str.join('_',[str(v) for v in variables]) - pylab.savefig(filename) # not good... - else: - pylab.show() - - return simul diff --git a/trash/dolo/numeric/newton_old.py b/trash/dolo/numeric/newton_old.py deleted file mode 100644 index 68411dd5..00000000 --- a/trash/dolo/numeric/newton_old.py +++ /dev/null @@ -1,177 +0,0 @@ -import numpy - -def simple_newton(f, x0, lb=None, ub=None, infos=False, verbose=False, maxit=50, tol=1e-8, eps=1e-8, numdiff=True): - '''Solves many independent systems f(x)=0 simultaneously using a simple gradient descent. - :param f: objective function to be solved with values p x N . The second output argument represents the derivative with - values in (p x p x N) - :param x0: initial value ( p x N ) - :return: solution x such that f(x) = 0 - ''' - - precision = x0.dtype # default tolerance should depend on precision - - from numpy.linalg import solve - - err = 1 - - it = 0 - while err > tol and it <= maxit: - - if not numdiff: - [res,dres] = f(x0) - else: - res = f(x0) - dres = numpy.zeros( (res.shape[0], x0.shape[0]), dtype=precision ) - for i in range(x0.shape[0]): - xi = x0.copy() - xi[i] += eps - resi = f(xi) - dres[:,i] = (resi - res)/eps - - dx = - solve(dres,res) - - x = x0 + dx - - print('x0 : {}'.format(x0)) - err = abs(res).max() - print('iteration {} {}'.format(it, err)) - - x0 = x - it += 1 - - if not infos: - return x - else: - return [x, it] - -def newton_solver(f, x0, lb=None, ub=None, infos=False, verbose=False, maxit=50, tol=1e-8, eps=1e-5, numdiff=False): - '''Solves many independent systems f(x)=0 simultaneously using a simple gradient descent. - :param f: objective function to be solved with values p x N . The second output argument represents the derivative with - values in (p x p x N) - :param x0: initial value ( p x N ) - :return: solution x such that f(x) = 0 - ''' - - precision = x0.dtype # default tolerance should depend on precision - - from dolo.numeric.serial_operations import serial_multiplication as stv, serial_solve - err = 1 - - it = 0 - while err > tol and it <= maxit: - if not numdiff: - [res,dres] = f(x0) - else: - res = f(x0) - dres = numpy.zeros( (res.shape[0], x0.shape[0], x0.shape[1]), dtype=precision ) - for i in range(x0.shape[0]): - xi = x0.copy() - xi[i,:] += eps - resi = f(xi) - dres[:,i,:] = (resi - res)/eps - - try: - dx = - serial_solve(dres,res) - except: - dx = - serial_solve(dres,res, debug=True) - x = x0 + dx - - err = abs(res).max() - - x0 = x - it += 1 - - if not infos: - return x - else: - return [x, it] - -def newton_solver_comp(f, x0, lb, ub, infos=False, backsteps=10, maxit=50, numdiff=False): - '''Solves many independent systems f(x)=0 simultaneously using a simple gradient descent. - :param f: objective function to be solved with values p x N . The second output argument represents the derivative with - values in (p x p x N) - :param x0: initial value ( p x N ) - :param lb: bounds for first variable - :param ub: bounds for second variable - :return: solution x such that f(x) = 0 - ''' - - from numpy import row_stack - - ind = x0.shape[0] - 1 - - def fun_lc(xx): - x = row_stack([xx, lb]) - res = f(x) - return res[:ind,:] - - def fun_uc(xx): - x = row_stack([xx, ub]) - res = f(x) - return res[:ind,:] - - [sol_nc, nit0] = newton_solver(f, x0, numdiff=True, infos=True) - lower_constrained = sol_nc[ind,:] < lb - upper_constrained = sol_nc[ind,:] > ub - not_constrained = - ( lower_constrained + upper_constrained ) - - - sol = sol_nc.copy() - -# sol[ind,:] = lb * lower_constrained + ub * upper_constrained + sol_nc[ind,:] * not_constrained -# nit = nit0 - - - [sol_lc, nit1] = newton_solver(fun_lc, x0[:-1,:], numdiff=True, infos=True) - [sol_uc, nit2] = newton_solver(fun_uc, x0[:-1,:], numdiff=True, infos=True) -# - nit = nit0 + nit1 + nit2 -# - sol_lc = row_stack([sol_lc, lb]) - sol_uc = row_stack([sol_uc, ub]) -# - lower_constrained = sol_nc[-1,:] < lb - upper_constrained = sol_nc[-1,:] > ub - not_constrained = - ( lower_constrained + upper_constrained ) -# - sol = sol_lc * lower_constrained + sol_uc * upper_constrained + sol_nc * not_constrained - - return [sol,nit] - - - - -from dolo.numeric.serial_operations import serial_inversion - -if __name__ == '__main__': - - p = 5 - N = 500 - - - import numpy.random - V = numpy.random.multivariate_normal([0]*p,numpy.eye(p),size=p) - print(V) - - M = numpy.zeros((p,p,N)) - for i in range(N): - M[:,:,i] = V - - from dolo.numeric.serial_operations import serial_multiplication as stm - - - MM = numpy.zeros( (p,N) ) - - - - import time - t = time.time() - for i in range(100): - T = serial_inversion(M) - s = time.time() - print('Elapsed :' + str(s-t)) - - - tt = stm(M,T) - for i in range(10): - print(tt[:,:,i]) diff --git a/trash/dolo/numeric/pea.py b/trash/dolo/numeric/pea.py deleted file mode 100644 index 90785aec..00000000 --- a/trash/dolo/numeric/pea.py +++ /dev/null @@ -1,156 +0,0 @@ -from numpy import tile, repeat -import numpy as np - -# compiled model is assumed to be a fgh model - -def compute_expectations_phi( cm, dr_x, s, x, p, nodes, weights): - - N = s.shape[1] - n_s = s.shape[0] - n_x = x.shape[0] - n_h = len(cm.model['variables_groups']['expectations']) - Q = len(weights) - - ee = nodes[numpy.newaxis,:,:] - - z = numpy.zeros( (n_h, N) ) - for i in range(Q): - e = ee[:,:,i:i+1] - w = weights[i] - S = cm.g(s,x,e,p) - X = dr_x(S) - h = cm.h(S,X,p) - z += h*w - return z - -def compute_expectations_psi( cm, dr_x, dr_z, s, x, p, nodes, weights): - - N = s.shape[1] - n_s = s.shape[0] - n_x = x.shape[0] - n_h = len(cm.model['variables_groups']['expectations']) - Q = len(weights) - - ss = tile(s, (1,Q)) - xx = tile(x, (1,Q)) - ee = repeat(nodes, N , axis=1) - SS = cm.g(ss,xx,ee,p) - - xxstart = dr_x(SS) - - ZZ = dr_z(SS) - - - ff = lambda xt: cm.f(SS,xt,ZZ,p) - [XX,nit] = newton_solver(ff , xxstart, numdiff=True,infos=True) - - hh = cm.h(SS,XX,p) - z = np.zeros( (n_x,N) ) - for i in range(Q): - z += weights[i] * hh[:,N*i:N*(i+1)] - - return z - - -from dolo.numeric.newton import newton_solver - - -def pea_solve( cm, grid, dr_x, dr_z, p, nodes, weights ): - - tol = 1e-9 - err = 10 - maxit = 500 - it = 0 - - x0 = dr_x(grid) - z0 = dr_z(grid) -# z0 = compute_expectations_phi( cm, dr_x, grid, x0, p, nodes, weights) - - while err > tol and it < maxit: - - it += 1 - - fobj = lambda x: cm.f( grid, x, z0, p ) - - x1 = newton_solver( fobj, x0, numdiff=True) - dr_x.set_values(x1) # I don't really need it - - z1 = compute_expectations_psi(cm, dr_x, dr_h, grid, x1, params, nodes, weights) - dr_z.set_values(z1) - - err = abs(z1 - z0).max() - err2 = abs(x1 - x0).max() - print(err, err2) - x0 = x1 - z0 = z1 - - print('finished in {} iterations'.format(it)) - - return [x0,z0] - - - -from dolo import yaml_import - -from trash.dolo.numeric.perturbations_to_states import approximate_controls - -model = yaml_import('../../../examples/models/compat/rbc_fgah.yaml') - - -dr_pert = approximate_controls(model, order=1, substitute_auxiliary=True) - -print(dr_pert.X_bar) - - - - -from dolo.numeric.interpolation.smolyak import SmolyakGrid -from dolo.numeric.time_iteration import time_iteration -dr_smol = time_iteration(model, smolyak_order=2, maxit=2, polish=True) - -smin = dr_smol.bounds[0,:] -smax = dr_smol.bounds[1,:] - - -dr_x = SmolyakGrid( smin, smax, 3) -dr_h = SmolyakGrid( smin, smax, 3) - -grid = dr_x.grid - -xh_init = dr_pert(dr_x.grid) - -n_x = len(model['variables_groups']['controls']) - - -dr_x.set_values( xh_init[:n_x,:] ) -dr_h.set_values( xh_init[n_x:,:] ) - - -import numpy - - -from dolo.compiler.cmodel_fgh import CModel_fgah, CModel_fgh_from_fgah - -cmt = CModel_fgah(model) - -cm = CModel_fgh_from_fgah(cmt) - - -sigma = cm.model.read_covariances() -params = cm.model.read_calibration()[2] -from dolo.numeric.quadrature import gauss_hermite_nodes -[nodes, weights] = gauss_hermite_nodes( [10], sigma) - -[x_sol, z_sol] = pea_solve( cm, grid, dr_x, dr_h, params, nodes, weights ) - -# - -dr_x.set_values(x_sol) - -# - - -# - -#dr_glob_3 = time_iteration(model_bis, smolyak_order=3) -#dr_glob_4 = time_iteration(model_bis, smolyak_order=4) diff --git a/trash/dolo/numeric/perturbations.py b/trash/dolo/numeric/perturbations.py deleted file mode 100644 index 0ceda895..00000000 --- a/trash/dolo/numeric/perturbations.py +++ /dev/null @@ -1,806 +0,0 @@ -from dolo.numeric.matrix_equations import second_order_solver, solve_sylvester -from dolo.numeric.tensor import sdot,mdot -from dolo.numeric.decision_rules import DynareDecisionRule as DDR -from dolo.numeric.decision_rules import TaylorExpansion - - -import numpy as np - -TOL = 1E-10 -TOL_RES = 1E-8 - -def solve_decision_rule(model, order=2, method='default',check_residuals=True, mlab=None, steady_state=None, solve_ss=False): - - - Sigma = model.get_distribution().sigma - y = np.array([model.calibration_dict[v] for v in model.variables]) - x = model.calibration['shocks'] - parms = model.calibration['parameters'] - - if steady_state != None: - y0 = steady_state - else: - y0 = y - - if solve_ss == True: - y = model.solve_for_steady_state(y0) - - - from dolo.compiler.compiler_statefree import PythonCompiler - pc = PythonCompiler(model) - - # pc = model.compiler - - if method == 'default': - p_dynamic = pc.compute_dynamic_pfile_cached(order,False,False) - elif method in ('sigma2','full'): - p_dynamic = pc.compute_dynamic_pfile_cached(order,False,True) - #p_dynamic = pc.compute_dynamic_pfile_cached(order,True,False) - else: - raise Exception('Unknown method : ' + method) - - # Build the approximation point to take the derivatives - yy = np.concatenate([y,y,y]) - xx = np.array(x) - xx = np.atleast_2d(xx) - parms = np.array(parms) - - derivatives = p_dynamic(yy,xx,parms) - - res = derivatives[0] - if check_residuals: - if abs(res).max() > TOL_RES: - print('Residuals\n') - for i,eq in enumerate(model.equations): - print('{0}\t:{1}\t:\t{2}'.format(i,res[i],eq)) - raise Exception("Initial values don't satisfy static equations") - - derivatives_ss = None - - if (method == 'sigma2'): - n_s = len(model.shocks) - n_v = len(model.variables) - nt = n_v*3 +n_s - nn = nt+1 - - if order == 2: - [f_0,f_1,f_2] = derivatives - ft_1 = f_1[:,:nt] - ft_2 = f_2[:,:nt,:nt] - f_ss = f_2[:,nt,nt] # derivatives with sigma^2 - derivatives = [f_0,ft_1,ft_2] - derivatives_ss = [f_ss] - - - elif order == 3: - [f_0,f_1,f_2,f_3] = derivatives - ft_1 = f_1[:,:nt] - ft_2 = f_2[:,:nt,:nt] - ft_3 = f_3[:,:nt,:nt,:nt] - f_ss = f_2[:,nt,nt] # derivatives with sigma^2 - f_1ss = f_3[:,:nt,nt,nt] - derivatives = [f_0,ft_1,ft_2,ft_3] - derivatives_ss = [f_ss,f_1ss] - - if method in ('sigma2','default'): - d = perturb_solver(derivatives, Sigma, max_order=order, derivatives_ss=derivatives_ss,mlab=mlab) - - elif method == 'full': - n_s = len(model.shocks) - n_v = len(model.variables) - n_p = len(model.parameters) - d = new_solver_with_p(derivatives,(n_v,n_s,n_p),max_order=order) - return d - - d['ys'] = np.array( y ) - d['Sigma'] = Sigma - -# return(TaylorExpansion(d)) - - return DDR( TaylorExpansion(d), model ) - - - -def perturb_solver(derivatives, Sigma, max_order=2, derivatives_ss=None, mlab=None): - - - if max_order == 1: - [f_0,f_1] = derivatives - elif max_order == 2: - [f_0,f_1,f_2] = derivatives - elif max_order == 3: - [f_0,f_1,f_2,f_3] = derivatives - else: - raise Exception('Perturbations not implemented at order {0}'.format(max_order)) - diff = derivatives - - f = diff - n = f[0].shape[0] # number of variables - s = f[1].shape[1] - 3*n - [n_v,n_s] = [n,s] - - - f1_A = f[1][:,:n] - f1_B = f[1][:,n:(2*n)] - f1_C = f[1][:,(2*n):(3*n)] - f1_D = f[1][:,(3*n):] - - ## first order - [ev,g_x] = second_order_solver(f1_A,f1_B,f1_C) - - res = np.dot(f1_A,np.dot(g_x,g_x)) + np.dot(f1_B,g_x) + f1_C - - mm = np.dot(f1_A, g_x) + f1_B - - g_u = - np.linalg.solve( mm , f1_D ) - - if max_order == 1: - d = {'ev':ev, 'g_a': g_x, 'g_e': g_u} - return d - - # we need it for higher order - V_a = np.concatenate( [np.dot(g_x,g_x),g_x,np.eye(n_v),np.zeros((s,n))] ) - V_e = np.concatenate( [np.dot(g_x,g_u),g_u,np.zeros((n_v,n_s)),np.eye(n_s)] ) - - # Translation - - f_1 = f[1] - f_2 = f[2] - f_d = f1_A - f_a = f1_B - g_a = g_x - g_e = g_u - n_a = n_v - n_e = n_s - - # Once for all ! - A = f_a + sdot(f_d,g_a) - B = f_d - C = g_a - A_inv = np.linalg.inv(A) - - - ################## - # Automatic code # - ################## - - #----------Computing order 2 - - order = 2 - - #--- Computing derivatives ('a', 'a') - - K_aa = + mdot(f_2,[V_a,V_a]) - L_aa = np.zeros( (n_v, n_v, n_v) ) - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aa + sdot(f_d,L_aa) - if mlab==None: - g_aa = solve_sylvester(A,B,C,D) - else: - n_d = D.ndim - 1 - n_v = C.shape[1] - CC = np.kron(C,C) - DD = D.reshape( n_v, n_v**n_d ) - [err,E] = mlab.gensylv(2,A,B,C,DD,nout=2) - g_aa = - E.reshape((n_v,n_v,n_v)) # check that - is correct - - - if order < max_order: - Y = L_aa + mdot(g_a,[g_aa]) + mdot(g_aa,[g_a,g_a]) - assert( abs(mdot(g_a,[g_aa]) - sdot(g_a,g_aa)).max() == 0) - Z = g_aa - V_aa = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('a', 'e') - - K_ae = + mdot(f_2,[V_a,V_e]) - L_ae = + mdot(g_aa,[g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ae) + K_ae - g_ae = - sdot(A_inv, const) - - if order < max_order: - Y = L_ae + mdot(g_a,[g_ae]) - Z = g_ae - V_ae = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('e', 'e') - - K_ee = + mdot(f_2,[V_e,V_e]) - L_ee = + mdot(g_aa,[g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ee) + K_ee - g_ee = - sdot(A_inv, const) - - if order < max_order: - Y = L_ee + mdot(g_a,[g_ee]) - Z = g_ee - V_ee = build_V(Y,Z,(n_a,n_e)) - - # manual - I = np.eye(n_v,n_v) - M_inv = np.linalg.inv( sdot(f1_A,g_a+I) + f1_B ) - K_ss = mdot(f_2[:,:n_v,:n_v],[g_e,g_e]) + sdot( f1_A, g_ee ) - rhs = - np.tensordot( K_ss, Sigma, axes=((1,2),(0,1)) ) #- mdot(h_2,[V_s,V_s]) - if derivatives_ss: - f_ss = derivatives_ss[0] - rhs -= f_ss - g_ss = sdot(M_inv,rhs) - ghs2 = g_ss/2 - - - - if max_order == 2: - d = { - 'ev': ev, - 'g_a': g_a, - 'g_e': g_e, - 'g_aa': g_aa, - 'g_ae': g_ae, - 'g_ee': g_ee, - 'g_ss': g_ss - } - return d - # /manual - - #----------Computing order 3 - - order = 3 - - #--- Computing derivatives ('a', 'a', 'a') - K_aaa = + 3*mdot(f_2,[V_a,V_aa]) + mdot(f_3,[V_a,V_a,V_a]) - L_aaa = + 3*mdot(g_aa,[g_a,g_aa]) - - #K_aaa = 2*( mdot(f_2,[V_aa,V_a]) ) + mdot(f_2,[V_a,V_aa]) + mdot(f_3,[V_a,V_a,V_a]) - #L_aaa = 2*( mdot(g_aa,[g_aa,g_a]) ) + mdot(g_aa,[g_a,g_aa]) - #K_aaa = ( mdot(f_2,[V_aa,V_a]) + mdot(f_2,[V_a,V_aa]) )*3.0/2.0 + mdot(f_3,[V_a,V_a,V_a]) - #L_aaa = ( mdot(g_aa,[g_aa,g_a]) + mdot(g_aa,[g_a,g_aa]) )*3.0/2.0 - - - #K_aaa = (K_aaa + K_aaa.swapaxes(3,2) + K_aaa.swapaxes(1,2) + K_aaa.swapaxes(1,2).swapaxes(2,3) + K_aaa.swapaxes(1,3) + K_aaa.swapaxes(1,3).swapaxes(2,3) )/6 - #L_aaa = (L_aaa + L_aaa.swapaxes(3,2) + L_aaa.swapaxes(1,2) + L_aaa.swapaxes(1,2).swapaxes(2,3) + L_aaa.swapaxes(1,3) + L_aaa.swapaxes(1,3).swapaxes(2,3) )/6 - - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aaa + sdot(f_d,L_aaa) - - - - if mlab == None: - g_aaa = solve_sylvester(A,B,C,D) - # this is much much faster - else: - n_d = D.ndim - 1 - n_v = C.shape[1] - CC = np.kron(np.kron(C,C),C) - DD = D.reshape( n_v, n_v**n_d ) - [err,E] = mlab.gensylv(3,A,B,C,DD,nout=2) - g_aaa = E.reshape((n_v,n_v,n_v,n_v)) - - #res = sdot(A,g_aaa) + sdot(B, mdot(g_aaa,[C,C,C])) - D - #print 'res : ' + str( abs(res).max() ) - - - if order < max_order: - Y = L_aaa + mdot(g_a,[g_aaa]) + mdot(g_aaa,[g_a,g_a,g_a]) - Z = g_aaa - V_aaa = build_V(Y,Z,(n_a,n_e)) - - # we transform g_aaa into a symmetric multilinear form - g_aaa = (g_aaa + g_aaa.swapaxes(3,2) + g_aaa.swapaxes(1,2) + g_aaa.swapaxes(1,2).swapaxes(2,3) + g_aaa.swapaxes(1,3) + g_aaa.swapaxes(1,3).swapaxes(2,3) )/6 - - #--- Computing derivatives ('a', 'a', 'e') - - K_aae = + mdot(f_2,[V_aa,V_e]) + 2*mdot(f_2,[V_a,V_ae]) + mdot(f_3,[V_a,V_a,V_e]) - L_aae = + mdot(g_aa,[g_aa,g_e]) + 2*mdot(g_aa,[g_a,g_ae]) + mdot(g_aaa,[g_a,g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aae) + K_aae - g_aae = - sdot(A_inv, const) - - if order < max_order: - Y = L_aae + mdot(g_a,[g_aae]) - Z = g_aae - V_aae = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('a', 'e', 'e') - - K_aee = + 2*mdot(f_2,[V_ae,V_e]) + mdot(f_2,[V_a,V_ee]) + mdot(f_3,[V_a,V_e,V_e]) - L_aee = + 2*mdot(g_aa,[g_ae,g_e]) + mdot(g_aa,[g_a,g_ee]) + mdot(g_aaa,[g_a,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aee) + K_aee - g_aee = - sdot(A_inv, const) - - if order < max_order: - Y = L_aee + mdot(g_a,[g_aee]) - Z = g_aee - V_aee = build_V(Y,Z,(n_a,n_e)) - - #--- Computing derivatives ('e', 'e', 'e') - - K_eee = + 3*mdot(f_2,[V_e,V_ee]) + mdot(f_3,[V_e,V_e,V_e]) - L_eee = + 3*mdot(g_aa,[g_e,g_ee]) + mdot(g_aaa,[g_e,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_eee) + K_eee - g_eee = - sdot(A_inv, const) - - if order < max_order: - Y = L_eee + mdot(g_a,[g_eee]) - Z = g_eee - V_eee = build_V(Y,Z,(n_a,n_e)) - - - #################################### - ## Compute sigma^2 correction term # - #################################### - - # ( a s s ) - - A = f_a + sdot(f_d,g_a) - I_e = np.eye(n_e) - - Y = g_e - Z = np.zeros((n_a,n_e)) - V_s = build_V(Y,Z,(n_a,n_e)) - - Y = mdot( g_ae, [g_a, I_e] ) - Z = np.zeros((n_a,n_a,n_e)) - V_as = build_V(Y,Z,(n_a,n_e)) - - Y = sdot(g_a,g_ss) + g_ss + np.tensordot(g_ee,Sigma) - Z = g_ss - V_ss = build_V(Y,Z,(n_a,n_e)) - - K_ass_1 = 2*mdot(f_2,[V_as,V_s] ) + mdot(f_3,[V_a,V_s,V_s]) - K_ass_1 = np.tensordot(K_ass_1,Sigma) - - K_ass_2 = mdot( f_2, [V_a,V_ss] ) - - K_ass = K_ass_1 + K_ass_2 - - L_ass = mdot( g_aa, [g_a, g_ss]) + np.tensordot( mdot(g_aee,[g_a, I_e, I_e]), Sigma) - - D = K_ass + sdot(f_d,L_ass) - - if derivatives_ss: - f_1ss = derivatives_ss[1] - D += mdot( f_1ss, [V_a]) + sdot( f1_A, mdot( g_aa, [ g_a , g_ss ]) ) - - g_ass = solve_sylvester( A, B, C, D) - - - # ( e s s ) - - A = f_a + sdot(f_d,g_a) - A_inv = np.linalg.inv(A) - I_e = np.eye(n_e) - - Y = g_e - Z = np.zeros((n_a,n_e)) - V_s = build_V(Y,Z,(n_a,n_e)) - - Y = mdot( g_ae, [g_e, I_e] ) - Z = np.zeros((n_a,n_e,n_e)) - V_es = build_V(Y,Z,(n_a,n_e)) - - Y = sdot(g_a,g_ss) + g_ss + np.tensordot(g_ee,Sigma) - Z = g_ss - V_ss = build_V(Y,Z,(n_a,n_e)) - - K_ess_1 = 2*mdot(f_2,[V_es,V_s] ) + mdot(f_3,[V_e,V_s,V_s]) - K_ess_1 = np.tensordot(K_ess_1,Sigma) - - K_ess_2 = mdot( f_2, [V_e,V_ss] ) - - K_ess = K_ess_1 + K_ess_2 - - L_ess = mdot( g_aa, [g_e, g_ss]) + np.tensordot( mdot(g_aee,[g_e, I_e, I_e]), Sigma) - L_ess += mdot( g_ass, [g_e]) - - D = K_ess + sdot(f_d,L_ess) - - g_ess = sdot( A_inv, -D) - - if max_order == 3: - d = {'ev':ev,'g_a':g_a,'g_e':g_e, 'g_aa':g_aa, 'g_ae':g_ae, 'g_ee':g_ee, - 'g_aaa':g_aaa, 'g_aae':g_aae, 'g_aee':g_aee, 'g_eee':g_eee, 'g_ss':g_ss, 'g_ass':g_ass,'g_ess':g_ess} - return d - - - -def new_solver_with_p(derivatives, sizes, max_order=2): - - if max_order == 1: - [f_0,f_1] = derivatives - elif max_order == 2: - [f_0,f_1,f_2] = derivatives - elif max_order == 3: - [f_0,f_1,f_2,f_3] = derivatives - diff = derivatives - - f = diff - #n = f[0].shape[0] # number of variables - #s = f[1].shape[1] - 3*n - [n_v,n_s,n_p] = sizes - - n = n_v - - f1_A = f[1][:,:n] - f1_B = f[1][:,n:(2*n)] - f1_C = f[1][:,(2*n):(3*n)] - f1_D = f[1][:,(3*n):((3*n)+n_s)] - f1_E = f[1][:,((3*n)+n_s):] - - ## first order - [ev,g_x] = second_order_solver(f1_A,f1_B,f1_C) - - mm = np.dot(f1_A, g_x) + f1_B - - g_u = - np.linalg.solve( mm , f1_D ) - g_p = - np.linalg.solve( mm , f1_E ) - - d = { - 'ev':ev, - 'g_a':g_x, - 'g_e':g_u, - 'g_p':g_p - } - - if max_order == 1: - return d - - # we need it for higher order - - V_x = np.concatenate( [np.dot(g_x,g_x),g_x,np.eye(n_v),np.zeros((n_s,n_v)), np.zeros((n_p,n_v))] ) - V_u = np.concatenate( [np.dot(g_x,g_u),g_u,np.zeros((n_v,n_s)),np.eye(n_s), np.zeros((n_p,n_s))] ) - V_p = np.concatenate( [np.dot(g_x,g_p),g_p,np.zeros((n_v,n_p)),np.zeros((n_s,n_p)), np.eye(n_p)] ) - V = [None, [V_x,V_u]] - - # Translation - n_a = n_v - n_e = n_s - n_p = g_p.shape[1] - - f_1 = f[1] - f_2 = f[2] - f_d = f1_A - f_a = f1_B - f_h = f1_C - f_u = f1_D - V_a = V_x - V_e = V_u - g_a = g_x - g_e = g_u - - # Once for all ! - A = f_a + sdot(f_d,g_a) - B = f_d - C = g_a - A_inv = np.linalg.inv(A) - - #----------Computing order 2 - - order = 2 - - #--- Computing derivatives ('a', 'a') - - K_aa = + mdot(f_2,[V_a,V_a]) - L_aa = np.zeros((n_v, n_a, n_a)) - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aa + sdot(f_d,L_aa) - g_aa = solve_sylvester(A,B,C,D) - - if order < max_order: - Y = L_aa + mdot(g_a,[g_aa]) + mdot(g_aa,[g_a,g_a]) - Z = g_aa - V_aa = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'e') - - K_ae = + mdot(f_2,[V_a,V_e]) - L_ae = + mdot(g_aa,[g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ae) + K_ae - g_ae = - sdot(A_inv, const) - - if order < max_order: - Y = L_ae + mdot(g_a,[g_ae]) - Z = g_ae - V_ae = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'p') - - K_ap = + mdot(f_2,[V_a,V_p]) - L_ap = + mdot(g_aa,[g_a,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ap) + K_ap - g_ap = - sdot(A_inv, const) - - if order < max_order: - Y = L_ap + mdot(g_a,[g_ap]) - Z = g_ap - V_ap = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'e') - - K_ee = + mdot(f_2,[V_e,V_e]) - L_ee = + mdot(g_aa,[g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ee) + K_ee - g_ee = - sdot(A_inv, const) - - if order < max_order: - Y = L_ee + mdot(g_a,[g_ee]) - Z = g_ee - V_ee = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'p') - - K_ep = + mdot(f_2,[V_e,V_p]) - L_ep = + mdot(g_aa,[g_e,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ep) + K_ep - g_ep = - sdot(A_inv, const) - - if order < max_order: - Y = L_ep + mdot(g_a,[g_ep]) - Z = g_ep - V_ep = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('p', 'p') - - K_pp = + mdot(f_2,[V_p,V_p]) - L_pp = + mdot(g_aa,[g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_pp) + K_pp - g_pp = - sdot(A_inv, const) - - if order < max_order: - Y = L_pp + mdot(g_a,[g_pp]) - Z = g_pp - V_pp = build_V(Y,Z,(n_a,n_e,n_p)) - - - d.update({ - 'g_aa':g_aa, - 'g_ae':g_ae, - 'g_ee':g_ee, - 'g_ap':g_ap, - 'g_ep':g_ep, - 'g_pp':g_pp - }) - if max_order == 2: - return d - - #----------Computing order 3 - - order = 3 - - #--- Computing derivatives ('a', 'a', 'a') - - K_aaa = + 3*mdot(f_2,[V_a,V_aa]) + mdot(f_3,[V_a,V_a,V_a]) - L_aaa = + 3*mdot(g_aa,[g_a,g_aa]) - - #We need to solve the infamous sylvester equation - #A = f_a + sdot(f_d,g_a) - #B = f_d - #C = g_a - D = K_aaa + sdot(f_d,L_aaa) - g_aaa = solve_sylvester(A,B,C,D) - - if order < max_order: - Y = L_aaa + mdot(g_a,[g_aaa]) + mdot(g_aaa,[g_a,g_a,g_a]) - Z = g_aaa - V_aaa = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'a', 'e') - - K_aae = + mdot(f_2,[V_aa,V_e]) + 2*mdot(f_2,[V_a,V_ae]) + mdot(f_3,[V_a,V_a,V_e]) - L_aae = + mdot(g_aa,[g_aa,g_e]) + 2*mdot(g_aa,[g_a,g_ae]) + mdot(g_aaa,[g_a,g_a,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aae) + K_aae - g_aae = - sdot(A_inv, const) - - if order < max_order: - Y = L_aae + mdot(g_a,[g_aae]) - Z = g_aae - V_aae = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'a', 'p') - - K_aap = + mdot(f_2,[V_aa,V_p]) + 2*mdot(f_2,[V_a,V_ap]) + mdot(f_3,[V_a,V_a,V_p]) - L_aap = + mdot(g_aa,[g_aa,g_p]) + 2*mdot(g_aa,[g_a,g_ap]) + mdot(g_aaa,[g_a,g_a,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aap) + K_aap - g_aap = - sdot(A_inv, const) - - if order < max_order: - Y = L_aap + mdot(g_a,[g_aap]) - Z = g_aap - V_aap = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'e', 'e') - - K_aee = + 2*mdot(f_2,[V_ae,V_e]) + mdot(f_2,[V_a,V_ee]) + mdot(f_3,[V_a,V_e,V_e]) - L_aee = + 2*mdot(g_aa,[g_ae,g_e]) + mdot(g_aa,[g_a,g_ee]) + mdot(g_aaa,[g_a,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aee) + K_aee - g_aee = - sdot(A_inv, const) - - if order < max_order: - Y = L_aee + mdot(g_a,[g_aee]) - Z = g_aee - V_aee = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'e', 'p') - ll = [ mdot(f_2,[V_ae,V_p]) , mdot(f_2,[V_ap,V_e]), mdot(f_2,[V_a,V_ep]) , mdot(f_3,[V_a,V_e,V_p]) ] - l = [ mdot(f_2,[V_ae,V_p]) , mdot(f_2,[V_ap,V_e]).swapaxes(2,3) , mdot(f_2,[V_a,V_ep]) , mdot(f_3,[V_a,V_e,V_p]) ] - - K_aep = + mdot(f_2,[V_ae,V_p]) + mdot(f_2,[V_ap,V_e]).swapaxes(2,3) + mdot(f_2,[V_a,V_ep]) + mdot(f_3,[V_a,V_e,V_p]) - L_aep = + mdot(g_aa,[g_ae,g_p]) + mdot(g_aa,[g_ap,g_e]).swapaxes(2,3) + mdot(g_aa,[g_a,g_ep]) + mdot(g_aaa,[g_a,g_e,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_aep) + K_aep - g_aep = - sdot(A_inv, const) - - if order < max_order: - Y = L_aep + mdot(g_a,[g_aep]) - Z = g_aep - V_aep = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('a', 'p', 'p') - - K_app = + 2*mdot(f_2,[V_ap,V_p]) + mdot(f_2,[V_a,V_pp]) + mdot(f_3,[V_a,V_p,V_p]) - L_app = + 2*mdot(g_aa,[g_ap,g_p]) + mdot(g_aa,[g_a,g_pp]) + mdot(g_aaa,[g_a,g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_app) + K_app - g_app = - sdot(A_inv, const) - - if order < max_order: - Y = L_app + mdot(g_a,[g_app]) - Z = g_app - V_app = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'e', 'e') - - K_eee = + 3*mdot(f_2,[V_e,V_ee]) + mdot(f_3,[V_e,V_e,V_e]) - L_eee = + 3*mdot(g_aa,[g_e,g_ee]) + mdot(g_aaa,[g_e,g_e,g_e]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_eee) + K_eee - g_eee = - sdot(A_inv, const) - - if order < max_order: - Y = L_eee + mdot(g_a,[g_eee]) - Z = g_eee - V_eee = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'e', 'p') - - K_eep = + mdot(f_2,[V_ee,V_p]) + 2*mdot(f_2,[V_e,V_ep]) + mdot(f_3,[V_e,V_e,V_p]) - L_eep = + mdot(g_aa,[g_ee,g_p]) + 2*mdot(g_aa,[g_e,g_ep]) + mdot(g_aaa,[g_e,g_e,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_eep) + K_eep - g_eep = - sdot(A_inv, const) - - if order < max_order: - Y = L_eep + mdot(g_a,[g_eep]) - Z = g_eep - V_eep = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('e', 'p', 'p') - - K_epp = + 2*mdot(f_2,[V_ep,V_p]) + mdot(f_2,[V_e,V_pp]) + mdot(f_3,[V_e,V_p,V_p]) - L_epp = + 2*mdot(g_aa,[g_ep,g_p]) + mdot(g_aa,[g_e,g_pp]) + mdot(g_aaa,[g_e,g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_epp) + K_epp - g_epp = - sdot(A_inv, const) - - if order < max_order: - Y = L_epp + mdot(g_a,[g_epp]) - Z = g_epp - V_epp = build_V(Y,Z,(n_a,n_e,n_p)) - - #--- Computing derivatives ('p', 'p', 'p') - - K_ppp = + 3*mdot(f_2,[V_p,V_pp]) + mdot(f_3,[V_p,V_p,V_p]) - L_ppp = + 3*mdot(g_aa,[g_p,g_pp]) + mdot(g_aaa,[g_p,g_p,g_p]) - - #We solve A*X + const = 0 - const = sdot(f_d,L_ppp) + K_ppp - g_ppp = - sdot(A_inv, const) - - if order < max_order: - Y = L_ppp + mdot(g_a,[g_ppp]) - Z = g_ppp - V_ppp = build_V(Y,Z,(n_a,n_e,n_p)) - - - d.update({ - 'g_aaa':g_aaa, - 'g_aae':g_aae, - 'g_aee':g_aee, - 'g_eee':g_eee, - 'g_aap':g_aap, - 'g_aep':g_aep, - 'g_eep':g_eep, - 'g_app':g_app, - 'g_epp':g_epp, - 'g_ppp':g_ppp - }) - - return d - - -def sigma2_correction(d,order,n_v,n_s,f1_A,f1_B, f_2, h_2, f_ss, Sigma): - - resp = dict() - - if order >= 2: - g_a = d['g_a'] - g_e = d['g_e'] - g_aa = d['g_aa'] - g_ae = d['g_ae'] - g_ee = d['g_ee'] - I = np.eye(n_v,n_v) - A_inv = np.linalg.inv( sdot(f1_A,g_a+I) + f1_B ) - V_s = np.zeros( (n_v*3+n_s+1,)) - V_s[-1] = 1 - K_ss = mdot(f_2[:,:n_v,:n_v],[g_e,g_e]) + sdot( f1_A, g_ee ) - rhs = np.tensordot( K_ss,Sigma) + mdot(h_2,[V_s,V_s]) - sigma2 = sdot(A_inv,-rhs) / 2 - - ds2 = sigma2 - d['g_ss'] - resp['sigma2'] = sigma2 - - if order == 3: - A = sdot(f1_A,g_a) + f1_B - A_inv = np.linalg.inv( A ) - B = f1_A - C = g_a - V_x = np.row_stack( [ - np.dot(g_a,g_a), - g_a, - I, - np.zeros((n_s,n_v)) - ]) - - D = mdot( f_ss, [V_x]) + sdot( f1_A, mdot( g_aa, [ g_a , ds2 ]) ) - dsigma2 = solve_sylvester(A,B,C,D) - resp['dsigma2'] = dsigma2 - resp['D'] = D - resp['f_ss'] = f_ss - return resp - -def build_V(X,Y,others): - - ddim = X.shape[1:] - new_dims = [(s,) + ddim for s in others] - l = [X,Y] + [np.zeros(d) for d in new_dims] - return np.concatenate(l, axis=0) diff --git a/trash/dolo/numeric/perturbations_dynare.py b/trash/dolo/numeric/perturbations_dynare.py deleted file mode 100644 index 2a9e56b6..00000000 --- a/trash/dolo/numeric/perturbations_dynare.py +++ /dev/null @@ -1,31 +0,0 @@ -# Here we use Dynare to compute the perturbations - - - - - - -def retrieve_DDR_from_matlab(name,mlab): - mlab.execute( 'drn = reorder_dr({0});'.format(name) ) - rdr = retrieve_from_matlab('drn',mlab) - ys = rdr['ys'][0,0].flatten() - ghx = rdr['ghx'][0,0] - ghu = rdr['ghu'][0,0] - [n_v,n_states] = ghx.shape - n_shocks = ghu.shape[1] - ghxx = rdr['ghxx'][0,0].reshape( (n_v,n_states,n_states) ) - ghxu = rdr['ghxu'][0,0].reshape( (n_v,n_states,n_shocks) ) - ghuu = rdr['ghuu'][0,0].reshape( (n_v,n_shocks,n_shocks) ) - ghs2 = rdr['ghs2'][0,0].flatten() - ddr = DDR( [ ys,[ghx,ghu],[ghxx,ghxu,ghuu] ] , ghs2 = ghs2 ) - return [ddr,rdr] - - -dr1_info_code = { -1: "the model doesn't define current variables uniquely.", -2: "problem in mjdgges.dll info(2) contains error code.", -3: "BK order condition not satisfied info(2) contains 'distance' absence of stable trajectory.", -4: "BK order condition not satisfied info(2) contains 'distance' indeterminacy.", -5: "BK rank condition not satisfied.", -6: "The jacobian matrix evaluated at the steady state is complex." -} diff --git a/trash/dolo/numeric/perturbations_to_states.py b/trash/dolo/numeric/perturbations_to_states.py deleted file mode 100644 index 6bc2c133..00000000 --- a/trash/dolo/numeric/perturbations_to_states.py +++ /dev/null @@ -1,367 +0,0 @@ -from dolo.numeric.decision_rules_states import CDR - -def approximate_controls(model, order=1, lambda_name=None, return_dr=True, verbose=True): - - from dolo.compiler.converter import model_to_fg - [g_fun, f_fun] = model_to_fg(model, order=order) - - # get steady_state - import numpy - - - - states = model.symbols_s['states'] - controls = model.symbols_s['controls'] - - parms = model.calibration['parameters'] - sigma = model.calibration['covariances'] - states_ss = model.calibration['states'] - controls_ss = model.calibration['controls'] - shocks_ss = model.calibration['shocks'] - - f_args_ss = numpy.concatenate( [states_ss, controls_ss, states_ss, controls_ss] ) - g_args_ss = numpy.concatenate( [states_ss, controls_ss, shocks_ss] ) - - f = f_fun( f_args_ss, parms) - g = g_fun( g_args_ss, parms) - - if lambda_name: - epsilon = 0.001 - sigma_index = [p.name for p in model.parameters].index(lambda_name) - pert_parms = parms.copy() - pert_parms[sigma_index] += epsilon - - g_pert = g_fun(g_args_ss, pert_parms) - sig2 = (g_pert[0] - g[0])/epsilon*2 - sig2_s = (g_pert[1] - g[1])/epsilon*2 - pert_sol = state_perturb(f, g, sigma, sigma2_correction = [sig2, sig2_s] ) - - else: - g = g_fun( g_args_ss, parms) - pert_sol = state_perturb(f, g, sigma ) - - n_s = len(states_ss) - n_c = len(controls_ss) - - if order == 1: - if return_dr: - S_bar = numpy.array( states_ss ) - X_bar = numpy.array( controls_ss ) - # add transitions of states to the d.r. - - - - X_s = pert_sol[0] - A = g[1][:,:n_s] + numpy.dot( g[1][:,n_s:n_s+n_c], X_s ) - B = g[1][:,n_s+n_c:] - dr = CDR([S_bar, X_bar, X_s]) - dr.A = A - dr.B = B - dr.sigma = sigma - return dr - return [controls_ss] + pert_sol - - if order == 2: - [[X_s,X_ss],[X_tt]] = pert_sol - X_bar = controls_ss + X_tt/2 - if return_dr: - S_bar = states_ss - S_bar = numpy.array(S_bar) - X_bar = numpy.array(X_bar) - dr = CDR([S_bar, X_bar, X_s, X_ss]) - A = g[1][:,:n_s] + numpy.dot( g[1][:,n_s:n_s+n_c], X_s ) - B = g[1][:,n_s+n_c:] - dr.sigma = sigma - dr.A = A - dr.B = B - return dr - return [X_bar, X_s, X_ss] - - - if order == 3: - [[X_s,X_ss,X_sss],[X_tt, X_stt]] = pert_sol - X_bar = controls_ss + X_tt/2 - X_s = X_s + X_stt/2 - if return_dr: - S_bar = states_ss - dr = CDR([S_bar, X_bar, X_s, X_ss, X_sss]) - dr.sigma = sigma - return dr - return [X_bar, X_s, X_ss, X_sss] - -def state_perturb(f_fun, g_fun, sigma, sigma2_correction=None, verbose=True): - """Computes a Taylor approximation of decision rules, given the supplied derivatives. - - The original system is assumed to be in the the form: - - .. math:: - - E_t f(s_t,x_t,s_{t+1},x_{t+1}) - - s_t = g(s_{t-1},x_{t-1}, \\lambda \\epsilon_t) - - where :math:`\\lambda` is a scalar scaling down the risk. the solution is a function :math:`\\varphi` such that: - - .. math:: - - x_t = \\varphi ( s_t, \\sigma ) - - The user supplies, a list of derivatives of f and g. - - :param f_fun: list of derivatives of f [order0, order1, order2, ...] - :param g_fun: list of derivatives of g [order0, order1, order2, ...] - :param sigma: covariance matrix of :math:`\\epsilon_t` - :param sigma2_correction: (optional) first and second derivatives of g w.r.t. sigma if :math:`g` explicitely depends - :math:`sigma` - - - Assuming :math:`s_t` , :math:`x_t` and :math:`\\epsilon_t` are vectors of size - :math:`n_s`, :math:`n_x` and :math:`n_x` respectively. - In general the derivative of order :math:`i` of :math:`f` is a multimensional array of size :math:`n_x \\times (N, ..., N)` - with :math:`N=2(n_s+n_x)` repeated :math:`i` times (possibly 0). - Similarly the derivative of order :math:`i` of :math:`g` is a multidimensional array of size :math:`n_s \\times (M, ..., M)` - with :math:`M=n_s+n_x+n_2` repeated :math:`i` times (possibly 0). - - - - """ - - import numpy as np - from numpy.linalg import solve - - approx_order = len(f_fun) - 1 # order of approximation - - [f0,f1] = f_fun[:2] - - [g0,g1] = g_fun[:2] - n_x = f1.shape[0] # number of controls - n_s = f1.shape[1]/2 - n_x # number of states - n_e = g1.shape[1] - n_x - n_s - n_v = n_s + n_x - - f_s = f1[:,:n_s] - f_x = f1[:,n_s:n_s+n_x] - f_snext = f1[:,n_v:n_v+n_s] - f_xnext = f1[:,n_v+n_s:] - - g_s = g1[:,:n_s] - g_x = g1[:,n_s:n_s+n_x] - g_e = g1[:,n_v:] - - A = np.row_stack([ - np.column_stack( [ np.eye(n_s), np.zeros((n_s,n_x)) ] ), - np.column_stack( [ -f_snext , -f_xnext ] ) - ]) - B = np.row_stack([ - np.column_stack( [ g_s, g_x ] ), - np.column_stack( [ f_s, f_x ] ) - ]) - - - - from dolo.numeric.extern.qz import qzordered - [S,T,Q,Z,eigval] = qzordered(A,B,n_s) - - # Check Blanchard=Kahn conditions - n_big_one = sum(eigval>1.0) - n_expected = n_x - if verbose: - print( "There are {} eigenvalues greater than 1. Expected: {}.".format( n_big_one, n_x ) ) - - if n_big_one != n_expected: - raise Exception("There should be exactly {} eigenvalues greater than one. Not {}.".format(n_x, n_big_one)) - - Q = Q.real # is it really necessary ? - Z = Z.real - - Z11 = Z[:n_s,:n_s] - Z12 = Z[:n_s,n_s:] - Z21 = Z[n_s:,:n_s] - Z22 = Z[n_s:,n_s:] - S11 = S[:n_s,:n_s] - T11 = T[:n_s,:n_s] - - # first order solution - C = solve(Z11.T, Z21.T).T - P = np.dot(solve(S11.T, Z11.T).T , solve(Z11.T, T11.T).T ) - Q = g_e - - if False: - from numpy import dot - test = f_s + dot(f_x,C) + dot( f_snext, g_s + dot(g_x,C) ) + dot(f_xnext, dot( C, g_s + dot(g_x,C) ) ) - print('Error: ' + str(abs(test).max())) - - if approx_order == 1: - return [C] - - # second order solution - from dolo.numeric.tensor import sdot, mdot - from numpy import dot - from dolo.numeric.matrix_equations import solve_sylvester - - f2 = f_fun[2] - g2 = g_fun[2] - g_ss = g2[:,:n_s,:n_s] - g_sx = g2[:,:n_s,n_s:n_v] - g_xx = g2[:,n_s:n_v,n_s:n_v] - - X_s = C - - - - V1_3 = g_s + dot(g_x,X_s) - V1 = np.row_stack([ - np.eye(n_s), - X_s, - V1_3, - dot( X_s, V1_3 ) - ]) - - K2 = g_ss + 2 * sdot(g_sx,X_s) + mdot(g_xx,[X_s,X_s]) - #L2 = - A = f_x + dot( f_snext + dot(f_xnext,X_s), g_x) - B = f_xnext - C = V1_3 - D = mdot(f2,[V1,V1]) + sdot(f_snext + dot(f_xnext,X_s),K2) - - X_ss = solve_sylvester(A,B,C,D) - -# test = sdot( A, X_ss ) + sdot( B, mdot(X_ss,[V1_3,V1_3]) ) + D - - - if not sigma == None: - g_ee = g2[:,n_v:,n_v:] - - v = np.row_stack([ - g_e, - dot(X_s,g_e) - ]) - - K_tt = mdot( f2[:,n_v:,n_v:], [v,v] ) - K_tt += sdot( f_snext + dot(f_xnext,X_s) , g_ee ) - K_tt += mdot( sdot( f_xnext, X_ss), [g_e, g_e] ) - K_tt = np.tensordot( K_tt, sigma, axes=((1,2),(0,1))) - - if sigma2_correction is not None: - K_tt += sdot( f_snext + dot(f_xnext,X_s) , sigma2_correction[0] ) - - L_tt = f_x + dot(f_snext, g_x) + dot(f_xnext, dot(X_s, g_x) + np.eye(n_x) ) - X_tt = solve( L_tt, - K_tt) - - if approx_order == 2: - if sigma == None: - return [X_s,X_ss] # here, we don't approximate the law of motion of the states - else: - return [[X_s,X_ss],[X_tt]] # here, we don't approximate the law of motion of the states - - # third order solution - - f3 = f_fun[3] - g3 = g_fun[3] - g_sss = g3[:,:n_s,:n_s,:n_s] - g_ssx = g3[:,:n_s,:n_s,n_s:n_v] - g_sxx = g3[:,:n_s,n_s:n_v,n_s:n_v] - g_xxx = g3[:,n_s:n_v,n_s:n_v,n_s:n_v] - - V2_3 = K2 + sdot(g_x,X_ss) - V2 = np.row_stack([ - np.zeros( (n_s,n_s,n_s) ), - X_ss, - V2_3, - dot( X_s, V2_3 ) + mdot(X_ss,[V1_3,V1_3]) - ]) - - K3 = g_sss + 3*sdot(g_ssx,X_s) + 3*mdot(g_sxx,[X_s,X_s]) + 2*sdot(g_sx,X_ss) - K3 += 3*mdot( g_xx,[X_ss,X_s] ) + mdot(g_xxx,[X_s,X_s,X_s]) - L3 = 3*mdot(X_ss,[V1_3,V2_3]) - - # A = f_x + dot( f_snext + dot(f_xnext,X_s), g_x) # same as before - # B = f_xnext # same - # C = V1_3 # same - D = mdot(f3,[V1,V1,V1]) + 3*mdot(f2,[ V2,V1 ]) + sdot(f_snext + dot(f_xnext,X_s),K3) - D += sdot( f_xnext, L3 ) - - X_sss = solve_sylvester(A,B,C,D) - - # now doing sigma correction with sigma replaced by l in the subscripts - - if not sigma is None: - g_se= g2[:,:n_s,n_v:] - g_xe= g2[:,n_s:n_v,n_v:] - - g_see= g3[:,:n_s,n_v:,n_v:] - g_xee= g3[:,n_s:n_v,n_v:,n_v:] - - - W_l = np.row_stack([ - g_e, - dot(X_s,g_e) - ]) - - I_e = np.eye(n_e) - - V_sl = g_se + mdot( g_xe, [X_s, np.eye(n_e)]) - - W_sl = np.row_stack([ - V_sl, - mdot( X_ss, [ V1_3, g_e ] ) + sdot( X_s, V_sl) - ]) - - K_ee = mdot(f3[:,:,n_v:,n_v:], [V1, W_l, W_l ]) - K_ee += 2 * mdot( f2[:,n_v:,n_v:], [W_sl, W_l]) - - # stochastic part of W_ll - - SW_ll = np.row_stack([ - g_ee, - mdot(X_ss, [g_e, g_e]) + sdot(X_s, g_ee) - ]) - - DW_ll = np.concatenate([ - X_tt, - dot(g_x, X_tt), - dot(X_s, sdot(g_x,X_tt )) + X_tt - ]) - - K_ee += mdot( f2[:,:,n_v:], [V1, SW_ll]) - - K_ = np.tensordot(K_ee, sigma, axes=((2,3),(0,1))) - - K_ += mdot(f2[:,:,n_s:], [V1, DW_ll]) - - def E(vec): - n = len(vec.shape) - return np.tensordot(vec,sigma,axes=((n-2,n-1),(0,1))) - - L = sdot(g_sx,X_tt) + mdot(g_xx,[X_s,X_tt]) - - L += E(g_see + mdot(g_xee,[X_s,I_e,I_e]) ) - - M = E( mdot(X_sss,[V1_3, g_e, g_e]) + 2*mdot(X_ss, [V_sl,g_e]) ) - M += mdot( X_ss, [V1_3, E( g_ee ) + sdot(g_x, X_tt)] ) - - - A = f_x + dot( f_snext + dot(f_xnext,X_s), g_x) # same as before - B = f_xnext # same - C = V1_3 # same - D = K_ + dot( f_snext + dot(f_xnext,X_s), L) + dot( f_xnext, M ) - - if sigma2_correction is not None: - g_sl = sigma2_correction[1][:,:n_s] - g_xl = sigma2_correction[1][:,n_s:(n_s+n_x)] - D += dot( f_snext + dot(f_xnext,X_s), g_sl + dot(g_xl,X_s) ) # constant - - X_stt = solve_sylvester(A,B,C,D) - - if approx_order == 3: - if sigma is None: - return [X_s,X_ss,X_sss] - else: - return [[X_s,X_ss,X_sss],[X_tt, X_stt]] - - -if __name__ == '__main__': - from dolo import * - model = yaml_import('examples/models/compat/rbc.yaml') - dr = approximate_controls(model) - print(dr.X_s) diff --git a/trash/dolo/numeric/portfolio_perturbation.py b/trash/dolo/numeric/portfolio_perturbation.py deleted file mode 100644 index 0e5d578c..00000000 --- a/trash/dolo/numeric/portfolio_perturbation.py +++ /dev/null @@ -1,208 +0,0 @@ -from trash.dolo.numeric import solver - - -def portfolios_to_deterministic(model,pf_names): - - ####### - ####### - - import re - regex = re.compile('.*<=(.*)<=.*') - for i,eq in enumerate(model.equations_groups['arbitrage']): - from trash.dolo.symbolic.symbolic import Variable, Equation - m = regex.match(eq.tags['complementarity']) - vs = m.group(1).strip() - if vs in pf_names: - v = Variable(vs) - neq = Equation(v,0) - neq.tag(**eq.tags) - model.equations_groups['arbitrage'][i] = neq - - print('Warning : initial model changed') - model.update() - - return model - -def solve_portfolio_model(model, pf_names, order=2, lambda_name='lam', guess=None): - - from dolo.compiler.compiler_python import GModel - if isinstance(model, GModel): - model = model.model - - pf_model = model - - from dolo import Variable, Parameter, Equation - import re - - n_states = len(pf_model.symbols_s['states']) - states = pf_model.symbols_s['states'] - steady_states = [Parameter(v.name+'_bar') for v in pf_model.symbols_s['states']] - n_pfs = len(pf_names) - - pf_vars = [Variable(v) for v in pf_names] - res_vars = [Variable('res_'+str(i)) for i in range(n_pfs)] - - - pf_parms = [Parameter('K_'+str(i)) for i in range(n_pfs)] - pf_dparms = [[Parameter('K_'+str(i)+'_'+str(j)) for j in range(n_states)] for i in range(n_pfs)] - - from sympy import Matrix - - # creation of the new model - - import copy - - new_model = copy.copy(pf_model) - - new_model.symbols_s['controls'] += res_vars - for v in res_vars + pf_vars: - new_model.calibration_s[v] = 0 - - - new_model.symbols_s['parameters'].extend(steady_states) - for p in pf_parms + Matrix(pf_dparms)[:]: - new_model.symbols_s['parameters'].append(p) - new_model.calibration_s[p] = 0 - - compregex = re.compile('(.*)<=(.*)<=(.*)') - - to_be_added_1 = [] - to_be_added_2 = [] - - expressions = Matrix(pf_parms) + Matrix(pf_dparms)*( Matrix(states) - Matrix(steady_states)) - - for n,eq in enumerate(new_model.equations_groups['arbitrage']): - if 'complementarity' in eq.tags: - tg = eq.tags['complementarity'] - [lhs,mhs,rhs] = compregex.match(tg).groups() - mhs = new_model.eval_string(mhs) - else: - mhs = None - if mhs in pf_vars: - i = pf_vars.index(mhs) - neq = Equation(mhs, expressions[i]) - neq.tag(**eq.tags) - eq_res = Equation(eq.gap, res_vars[i]) - eq_res.tag(eq_type='arbitrage') - to_be_added_2.append(eq_res) - new_model.equations_groups['arbitrage'][n] = neq - to_be_added_1.append(neq) - - # new_model.equations_groups['arbitrage'].extend(to_be_added_1) - new_model.equations_groups['arbitrage'].extend(to_be_added_2) - new_model.update() - - print("number of equations {}".format(len(new_model.equations))) - print("number of arbitrage equations {}".format( len(new_model.equations_groups['arbitrage'])) ) - - print('parameters_ordering') - print("number of parameters {}".format(new_model.symbols['parameters'])) - print("number of parameters {}".format(new_model.parameters)) - - - # now, we need to solve for the optimal portfolio coefficients - from trash.dolo.numeric.perturbations_to_states import approximate_controls - - dr = approximate_controls(new_model) - print('ok') - - import numpy - - n_controls = len(model.symbols_s['controls']) - - def constant_residuals(x, return_dr=False): - d = {} - for i in range(n_pfs): - p = pf_parms[i] - v = pf_vars[i] - d[p] = x[i] - d[v] = x[i] - new_model.set_calibration(d) - # new_model.parameters_values[p] = x[i] - # new_model.init_values[v] = x[i] - if return_dr: - dr = approximate_controls(new_model, order=1, return_dr=True, lambda_name='lam') - return dr - X_bar, X_s, X_ss = approximate_controls(new_model, order=2, return_dr=False, lambda_name="lam") - - return X_bar[n_controls-n_pfs:n_controls] - - - if guess is not None: - x0 = numpy.array(guess) - else: - x0 = numpy.zeros(n_pfs) - - print('Zero order portfolios') - print('Initial guess: {}'.format(x0)) - print('Initial error: {}'.format( constant_residuals(x0) )) - - portfolios_0 = solver(constant_residuals, x0) - print('Solution: {}'.format(portfolios_0)) - print('Final error: {}'.format( constant_residuals(portfolios_0) )) - - if order == 1: - dr = constant_residuals(portfolios_0, return_dr=True) - return dr - - def dynamic_residuals(X, return_dr=False): - x = X[:,0] - dx = X[:,1:] - d = {} - for i in range(n_pfs): - p = pf_parms[i] - v = pf_vars[i] - d[p] = x[i] - d[v] = x[i] - for j in range(n_states): - d[pf_dparms[i][j]] = dx[i,j] - new_model.set_calibration(d) - if return_dr: - dr = approximate_controls(new_model, order=2, lambda_name='lam') - return dr - else: - [X_bar, X_s, X_ss, X_sss] = approximate_controls(new_model, order=3, return_dr=False, lambda_name='lam') - crit = numpy.column_stack([ - X_bar[n_controls-n_pfs:n_controls], - X_s[n_controls-n_pfs:n_controls,:], - ]) - return crit - - - - y0 = numpy.column_stack([x0, numpy.zeros((n_pfs, n_states))]) - print('Initial error:') - err = (dynamic_residuals(y0)) - - print( abs(err).max() ) - portfolios_1 = solver(dynamic_residuals, y0) - - print('First order portfolios : ') - print(portfolios_1) - - print('Final error:') - print(dynamic_residuals(portfolios_1)) - - dr = dynamic_residuals(portfolios_1, return_dr=True) - - # TODO: remove coefficients of criteria - - return dr - -if __name__ == '__main__': - - from dolo import * - #model = yaml_import('examples/models/compat/open_economy_with_pf_pert.yaml') -# model = yaml_import('/home/pablo/Documents/Research/Thesis/chapter_1/code/capital.yaml') - model = yaml_import('/home/pablo/Documents/Research/Thesis/chapter_1/code/capital_pert.yaml') - - print(model) - print(model.calibration['covariances']) - - - sol = solve_portfolio_model(model,['x_1','x_2'], guess=[-0.3,-0.5], order=1) - print( sol.order ) - - print("Function has returned.") - print(sol) - print(sol.X_s) diff --git a/trash/dolo/numeric/risky_ss.py b/trash/dolo/numeric/risky_ss.py deleted file mode 100644 index fa847995..00000000 --- a/trash/dolo/numeric/risky_ss.py +++ /dev/null @@ -1,276 +0,0 @@ -# warning : this code is EXPERIMENTAL -# don't use it in serious projects unless you know what you are doing - - -from __future__ import division - -import numpy - -from dolo.numeric.decision_rules_states import CDR -from trash.dolo.numeric import solver - -numpy.set_printoptions(precision=3, suppress=3) - -import numpy as np - - -def solve_model_around_risky_ss(model, verbose=False, return_dr=True, initial_sol=None): - - - - if initial_sol == None: - if model.model_spec == 'portfolios': - print('This is a portfolio model ! Converting to deterministic one.') - from dolo.algos.portfolio_perturbation import portfolios_to_deterministic - model = portfolios_to_deterministic(model,['x_1','x_2']) - model.check() - from trash.dolo.numeric.perturbations_to_states import approximate_controls - perturb_sol = approximate_controls(model, order = 1, return_dr=False) - [X_bar,X_s] = perturb_sol - else: - perturb_sol = initial_sol - [X_bar,X_s] = perturb_sol - - - # reduce X_s to the real controls (remove auxiliary variables - - X_s = X_s[:len(model.symbols['controls']),:] - X_bar = X_bar[:len(model.symbols['controls'])] - X_bar = numpy.array(X_bar) - - if abs(X_s.imag).max() < 1e-10: - X_s = X_s.real - else: - raise( Exception('Complex decision rule') ) - - print('Perturbation solution found') - - #print model.parameters - #exit() - - X_bar_1 = X_bar - X_s_1 = X_s - - #model = yaml_import(filename) - - #from dolo.symbolic.symbolic import Parameter - - [S_bar, X_bar, X_s, P] = solve_risky_ss(model.model, X_bar, X_s, verbose=verbose) - - if return_dr: - cdr = CDR([S_bar, X_bar, X_s]) -# cdr.P = P - return cdr - return [S_bar, X_bar, X_s, P] - - - -def solve_risky_ss(model, X_bar, X_s, verbose=False): - - import numpy - from dolo.compiler.function_compiler import compile_function - import time - from dolo.compiler.compiler_functions import simple_global_representation - - - parms = model.calibration['parameters'] - sigma = model.calibration['covariances'] - - sgm = simple_global_representation(model) - - states = sgm['states'] - controls = sgm['controls'] - shocks = sgm['shocks'] - parameters = sgm['parameters'] - f_eqs = sgm['f_eqs'] - g_eqs = sgm['g_eqs'] - - - g_args = [s(-1) for s in states] + [c(-1) for c in controls] + shocks - f_args = states + controls + [v(1) for v in states] + [v(1) for v in controls] - p_args = parameters - - - - g_fun = compile_function(g_eqs, g_args, p_args, 2) - f_fun = compile_function(f_eqs, f_args, p_args, 3) - - - epsilons_0 = np.zeros((sigma.shape[0])) - - from numpy import dot - from dolo.numeric.tensor import sdot,mdot - - def residuals(X, sigma, parms, g_fun, f_fun): - - import numpy - - dummy_x = X[0:1,0] - X_bar = X[1:,0] - S_bar = X[0,1:] - X_s = X[1:,1:] - - [n_x,n_s] = X_s.shape - - n_e = sigma.shape[0] - - xx = np.concatenate([S_bar, X_bar, epsilons_0]) - - [g_0, g_1, g_2] = g_fun(xx, parms) - [f_0,f_1,f_2,f_3] = f_fun( np.concatenate([S_bar, X_bar, S_bar, X_bar]), parms) - - res_g = g_0 - S_bar - - # g is a first order function - g_s = g_1[:,:n_s] - g_x = g_1[:,n_s:n_s+n_x] - g_e = g_1[:,n_s+n_x:] - g_se = g_2[:,:n_s,n_s+n_x:] - g_xe = g_2[:, n_s:n_s+n_x, n_s+n_x:] - - - # S(s,e) = g(s,x,e) - S_s = g_s + dot(g_x, X_s) - S_e = g_e - S_se = g_se + mdot(g_xe,[X_s, numpy.eye(n_e)]) - - - # V(s,e) = [ g(s,x,e) ; x( g(s,x,e) ) ] - V_s = np.row_stack([ - S_s, - dot( X_s, S_s ) - ]) # *** - - V_e = np.row_stack([ - S_e, - dot( X_s, S_e ) - ]) - - V_se = np.row_stack([ - S_se, - dot( X_s, S_se ) - ]) - - # v(s) = [s, x(s)] - v_s = np.row_stack([ - numpy.eye(n_s), - X_s - ]) - - - # W(s,e) = [xx(s,e); yy(s,e)] - W_s = np.row_stack([ - v_s, - V_s - ]) - - #return - - nn = n_s + n_x - f_v = f_1[:,:nn] - f_V = f_1[:,nn:] - f_1V = f_2[:,:,nn:] - f_VV = f_2[:,nn:,nn:] - f_1VV = f_3[:,:,nn:,nn:] - - - # E = lambda v: np.tensordot(v, sigma, axes=((2,3),(0,1)) ) # expectation operator - - F = f_0 + 0.5*np.tensordot( mdot(f_VV,[V_e,V_e]), sigma, axes=((1,2),(0,1)) ) - - F_s = sdot(f_1, W_s) - f_see = mdot(f_1VV, [W_s, V_e, V_e]) + 2*mdot(f_VV, [V_se, V_e]) - F_s += 0.5 * np.tensordot(f_see, sigma, axes=((2,3),(0,1)) ) # second order correction - - resp = np.row_stack([ - np.concatenate([dummy_x,res_g]), - np.column_stack([F,F_s]) - ]) - - return resp - - # S_bar = s_fun_init( numpy.atleast_2d(X_bar).T ,parms).flatten() - # S_bar = S_bar.flatten() - S_bar = model.calibration['states'] - S_bar = np.array(S_bar) - - X0 = np.row_stack([ - np.concatenate([np.zeros(1),S_bar]), - np.column_stack([X_bar,X_s]) - ]) - - fobj = lambda X: residuals(X, sigma, parms, g_fun, f_fun) - - if verbose: - val = fobj(X0) - print('val') - print(val) - - # exit() - - t = time.time() - - sol = solver(fobj,X0, method='lmmcp', verbose=verbose, options={'preprocess':False, 'eps1':1e-15, 'eps2': 1e-15}) - - if verbose: - print('initial guess') - print(X0) - print('solution') - print sol - print('initial residuals') - print(fobj(X0)) - print('residuals') - print fobj(sol) - s = time.time() - - - if verbose: - print('Elapsed : {0}'.format(s-t)) - #sol = solver(fobj,X0, method='fsolve', verbose=True, options={'preprocessor':False}) - - norm = lambda x: numpy.linalg.norm(x,numpy.inf) - if verbose: - print( "Initial error: {0}".format( norm( fobj(X0)) ) ) - print( "Final error: {0}".format( norm( fobj(sol) ) ) ) - - print("Solution") - print(sol) - - X_bar = sol[1:,0] - S_bar = sol[0,1:] - X_s = sol[1:,1:] - - # compute transitions - n_s = len(states) - n_x = len(controls) - [g, dg, junk] = g_fun( np.concatenate( [S_bar, X_bar, epsilons_0] ), parms) - g_s = dg[:,:n_s] - g_x = dg[:,n_s:n_s+n_x] - - P = g_s + dot(g_x, X_s) - - - if verbose: - eigenvalues = numpy.linalg.eigvals(P) - print eigenvalues - eigenvalues = [abs(e) for e in eigenvalues] - eigenvalues.sort() - print(eigenvalues) - - return [S_bar, X_bar, X_s, P] - - -# -# -# -#fname = 'capital' -# -#filename = 'models/{0}.yaml'.format(fname) -# -#[S_bar, X_bar, X_s, P] = solve_model_around_risky_ss(filename) -if __name__ == '__main__': - from dolo import yaml_import - fname = '/home/pablo/Documents/Research/Thesis/chapter_4/code/models/open_economy.yaml' - #fname = 'examples/models/compat/open_economy_with_pf.yaml' - model = yaml_import(fname) - sol = solve_model_around_risky_ss(model, verbose=True) \ No newline at end of file diff --git a/trash/dolo/numeric/serial_operations_cython.pyx b/trash/dolo/numeric/serial_operations_cython.pyx deleted file mode 100644 index daa814d7..00000000 --- a/trash/dolo/numeric/serial_operations_cython.pyx +++ /dev/null @@ -1,125 +0,0 @@ -import numpy as np -cimport numpy as np -from numpy.linalg import inv - - -DTYPE = np.float64 -ctypedef np.float64_t DTYPE_t - -cimport cython - - - -@cython.boundscheck(False) -@cython.wraparound(False) -def serial_multiplication(np.ndarray[DTYPE_t, ndim=3] A, np.ndarray[DTYPE_t, ndim=3] B): - - cdef int i,k,j,n,I,K,J,N - I = A.shape[0] - J = A.shape[1] - N = A.shape[2] - K = B.shape[1] - - cdef np.ndarray[DTYPE_t, ndim=3] C = np.zeros( (I,K,N), dtype = np.float64 ) - - cdef double* A_d = A.data - cdef double* B_d = B.data - cdef double* C_d = C.data - - cdef int rg, rm, rd - - for i in range(I): - for k in range(K): - for j in range(J): - rg = i*N*K + k*N - rm = i*N*J + j*N - rd = j*N*K + k*N - for n in range(N): - C_d[rg] += A_d[rm]*B_d[rd] - rg +=1 - rm +=1 - rd +=1 - - return C - - - - -@cython.boundscheck(False) -@cython.wraparound(False) -def serial_dot_21(np.ndarray[DTYPE_t, ndim=3] A, np.ndarray[DTYPE_t, ndim=2] B): - - cdef int i,k,j,n,I,K,J,N - I = A.shape[0] - J = A.shape[1] - N = A.shape[2] -# K = B.shape[1] : B.shape should be [J,N] - - # assert(B.shape[0]==J) - # assert(B.shape[2]==N) - - cdef np.ndarray[DTYPE_t, ndim=2] resp - # cdef np.ndarray[DTYPE_t, ndim=1] T - - resp = np.zeros( (I,N), dtype = np.float64 ) - for i in range(I): - for j in range(J): - for n in range(N): - resp[i,n] += A[i,j,n]*B[j,n] - - return resp - -@cython.boundscheck(False) -@cython.wraparound(False) -def serial_dot_11(np.ndarray[DTYPE_t, ndim=2] A, np.ndarray[DTYPE_t, ndim=2] B): - - cdef int i,k,j,n,I,K,J,N - I = A.shape[0] - N = A.shape[1] - # K = B.shape[1] : B.shape should be [I,N] - - assert(B.shape[0]==I) - assert(B.shape[1]==N) - - cdef np.ndarray[DTYPE_t, ndim=2] resp - # cdef np.ndarray[DTYPE_t, ndim=1] T - - resp = np.zeros( (I,N), dtype = np.float64 ) - for i in range(I): - for n in range(N): - resp[i,n] += A[i,n]*B[i,n] - - return resp - - -@cython.boundscheck(False) -@cython.wraparound(False) -def serial_inversion(np.ndarray[DTYPE_t, ndim=3] M): - ''' - - :param M: a pxpxN array - :return: a pxpxN array T such that T(:,:,i) is the inverse of M(:,:,i) - ''' - - import numpy - from numpy.linalg import inv - - cdef np.ndarray[DTYPE_t, ndim=3] T - cdef np.ndarray[DTYPE_t, ndim=2] tmp - cdef int i - - p = M.shape[0] - assert(M.shape[1] == p) - N = M.shape[2] - - M = np.asfortranarray(M) - - T = np.zeros((p,p,N),order='F') - - for i in range(N): - tmp = M[:,:,i] - T[:,:,i] = inv(tmp) - - T = numpy.ascontiguousarray(T) - - return T \ No newline at end of file diff --git a/trash/dolo/numeric/solver.py b/trash/dolo/numeric/solver.py deleted file mode 100644 index 9ae4bb35..00000000 --- a/trash/dolo/numeric/solver.py +++ /dev/null @@ -1,129 +0,0 @@ -import numpy as np - - -def solver(fobj, x0, lb=None, ub=None, jac=None, method='lmmcp', infos=False, serial_problem=False, verbose=False, options={}): - - - in_shape = x0.shape - - if serial_problem: - ffobj = fobj - else: - ffobj = lambda x: fobj(x.reshape(in_shape)).flatten() - - - # standardize jacobian - if jac is not None: - if not serial_problem: - pp = np.prod(in_shape) - def Dffobj(t): - tt = t.reshape(in_shape) - dval = jac(tt) - return dval.reshape( (pp,pp) ) - else: - Dffobj = jac - elif serial_problem: - from dolo.numeric.newton import SerialDifferentiableFunction - Dffobj = SerialDifferentiableFunction(fobj, in_shape) - else: - Dffobj = MyJacobian(ffobj) - - - if lb == None: - lb = -np.inf*np.ones(len(x0.flatten())) - if ub == None: - ub = np.inf*np.ones(len(x0.flatten())).flatten() - - if not serial_problem: - lb = lb.flatten() - ub = ub.flatten() - - if not serial_problem: - x0 = x0.flatten() - - - if method == 'fsolve': - import scipy.optimize as optimize - factor = options.get('factor') - factor = factor if factor else 1 - [sol,infodict,ier,msg] = optimize.fsolve(ffobj, x0, fprime=Dffobj, factor=factor, full_output=True, xtol=1e-10, epsfcn=1e-9) - if ier != 1: - print(msg) - - elif method == 'newton': - from dolo.numeric.newton import serial_newton as newton_solver - fun = lambda x: [ffobj(x), Dffobj(x) ] - [sol,nit] = newton_solver(fun,x0, verbose=verbose) - - elif method == 'lmmcp': - from dolo.numeric.extern.lmmcp import lmmcp - sol = lmmcp(lambda t: -ffobj(t), lambda u: -Dffobj(u),x0,lb,ub,verbose=verbose,options=options) - - elif method == 'ncpsolve': - from dolo.numeric.ncpsolve import ncpsolve - fun = lambda x: [ffobj(x), Dffobj(x) ] - if serial_problem: - jactype = 'serial' - [sol,nit] = ncpsolve(fun,lb,ub,x0, verbose=verbose, infos=True, jactype='serial') - - else: - raise Exception('Unknown method : '+str(method)) - sol = sol.reshape(in_shape) - - if infos: - return [sol, nit] - else: - return sol - - -def MyJacobian(fun, eps=1e-6): - - def rfun(x): - n = len(x) - x0 = x.copy() - y0 = fun(x) - D = np.zeros( (len(y0),len(x0)) ) - for i in range(n): - delta = np.zeros(len(x)) - delta[i] = eps - y1 = fun(x+delta) - y2 = fun(x-delta) - D[:,i] = (y1 - y2)/eps/2 - return D - return rfun - - - -def MySerialJacobian(fun, shape, eps=1e-6): - - def rfun(x): - - x = x.reshape(shape) - - # x0 = x.copy() - p = x.shape[1] - N = x.shape[0] - - y0 = fun(x) - - assert( y0.shape[1] == p) - assert( y0.shape[0] == N) - - Dc = np.zeros( (N,p,p) ) # compressed jacobian - for i in range(p): - delta = np.zeros((N,p)) - delta[:,i] = eps - y1 = fun(x+delta) - y2 = fun(x-delta) - Dc[:,:,i] = (y1 - y2)/eps/2 - - return Dc.swapaxes(0,1) - -# D = np.zeros((p,N,p,N)) -# for n in range(N): -# D[:,n,:,n] = Dc[:,:,n].T -# -# return D.reshape(p*N, p*N) -# #return D - - return rfun diff --git a/trash/dolo/tests/test_bigsystem.py b/trash/dolo/tests/test_bigsystem.py deleted file mode 100644 index 558bd606..00000000 --- a/trash/dolo/tests/test_bigsystem.py +++ /dev/null @@ -1,20 +0,0 @@ -def test_big_system(): - from dolo import yaml_import - from dolo.algos.dtcscc.nonlinearsystem import nonlinear_system - - from dolo.algos.dtcscc.time_iteration import time_iteration - - model = yaml_import('examples/models/compat/rbc.yaml') - - dr = time_iteration(model, grid={'type': 'smolyak', 'mu': 3}, verbose=True) - sol = nonlinear_system(model, grid={'type': 'smolyak', 'mu': 3}) - - diff = (sol.__values__) - dr.__values__ - assert(abs(diff).max() < 1e-6) - - sol_high_precision = nonlinear_system(model, grid={'type': 'smolyak', 'mu': 5}, dr0=sol) - - assert(sol_high_precision.grid.shape == (145, 2)) - -if __name__ == '__main__': - test_big_system() diff --git a/trash/dolo/tests/test_dtcscc_direct.py b/trash/dolo/tests/test_dtcscc_direct.py deleted file mode 100644 index e52ff846..00000000 --- a/trash/dolo/tests/test_dtcscc_direct.py +++ /dev/null @@ -1,31 +0,0 @@ -def test_direct(): - - from dolo import yaml_import - from dolo.algos.dtcscc.perturbations import approximate_controls - from dolo.algos.dtcscc.time_iteration import time_iteration_direct, time_iteration - - model = yaml_import("examples/models/compat/rbc_full.yaml") - - # Check without complementarity conditions - dr = time_iteration_direct(model, with_complementarities=False) - ddr = time_iteration(model, with_complementarities=False) - - x0 = dr(dr.grid) - x1 = ddr(dr.grid) - - print(abs(x1 - x0).max()<1e-5) - - - # Check with complementarity conditions - dr = time_iteration_direct(model, with_complementarities=True) - ddr = time_iteration(model, with_complementarities=True) - - x0 = dr(dr.grid) - x1 = ddr(dr.grid) - - print(abs(x1 - x0).max()<1e-5) - - -if __name__ =='__main__': - - test_direct() diff --git a/trash/dolo/tests/test_dynare_model.py b/trash/dolo/tests/test_dynare_model.py deleted file mode 100644 index 5efce9d9..00000000 --- a/trash/dolo/tests/test_dynare_model.py +++ /dev/null @@ -1,37 +0,0 @@ -from dolo import * - -from dolo.algos.dynare.perturbations import * -from dolo.algos.dynare.simulations import * - -def test_dynare_model(): - - model = yaml_import("examples/models/compat/rbc_dynare.yaml") - - dr = solve_decision_rule(model, order=2) - print(dr['ys']) - - n_a = len(model.variables) + len(model.symbols['shocks']) - - irf = simulate(dr, n_exp=0, start='risky') - print(irf) - assert(irf.shape == (41,n_a)) - assert(irf.std().max()<0.00001) - - irf = simulate(dr, n_exp=0, start='deterministic') - assert(irf.shape == (41,n_a)) - assert(irf.std().max()>0.00001) - - sim = simulate(dr, n_exp=1) - assert(sim.shape == (41,n_a)) - assert(sim.std().max()>0.00001) - - sim2 = simulate(dr, n_exp=10) - assert(sim2.shape == (10,41,n_a)) - - sim_shock = simulate(dr, n_exp=0, shock=[0.01]) - - - -if __name__ == '__main__': - - test_dynare_model() diff --git a/trash/dolo/tests/test_errors.py b/trash/dolo/tests/test_errors.py deleted file mode 100644 index 29f375cb..00000000 --- a/trash/dolo/tests/test_errors.py +++ /dev/null @@ -1,67 +0,0 @@ -def test_omega_errors(): - - from dolo import yaml_import - from dolo.algos.dtcscc.time_iteration import time_iteration as time_iteration - - model = yaml_import('examples/models/compat/rbc.yaml') - - from dolo.algos.dtcscc.perturbations import approximate_controls - - dr = approximate_controls(model) - dr_global = time_iteration(model, verbose=False, pert_order=1) - - sigma = model.get_distribution().sigma - - s_0 = dr.S_bar - - from dolo.algos.dtcscc.accuracy import omega - - res_1 = omega(model, dr, grid=dict(orders=[10, 10]), time_discount=0.96) - res_2 = omega(model, dr_global) - - print(res_1) - print(res_2) - - -def test_denhaan_errors(): - - from dolo import yaml_import - from dolo.algos.dtcscc.time_iteration import time_iteration as time_iteration - - model = yaml_import('examples/models/compat/rbc.yaml') - - from dolo.algos.dtcscc.perturbations import approximate_controls - - dr = approximate_controls(model) - - dr_global = time_iteration(model, verbose=False) - - s0 = model.calibration['states'] - sigma = model.get_distribution().sigma - - from dolo.algos.dtcscc.accuracy import denhaanerrors - - denerr_1 = denhaanerrors(model, dr) - denerr_2 = denhaanerrors(model, dr_global) - - print(denerr_1) - print(denerr_2) - print(denerr_2['max_errors'][0]) - - - assert( max(denerr_2['max_errors']) < 10-7) # errors with solyak colocations at order 4 are very small - # assert( max(denerr_1['mean_errors']) < 10-7) - -if __name__ == '__main__': - test_omega_errors() - test_denhaan_errors() - - - - - - -if __name__ == '__main__': - - test_denhaan_errors() - test_omega_errors() diff --git a/trash/dolo/tests/test_fg_higher_order_perturbations.py b/trash/dolo/tests/test_fg_higher_order_perturbations.py deleted file mode 100644 index da3fb583..00000000 --- a/trash/dolo/tests/test_fg_higher_order_perturbations.py +++ /dev/null @@ -1,17 +0,0 @@ -def test_fga_higher_order_perturbations(): - - from dolo import yaml_import - from dolo.algos.dtcscc.perturbations_higher_order import approximate_controls - - model = yaml_import('examples/models/compat/rbc.yaml') - # for i in [1,2,3]: - dr1 = approximate_controls(model, order=1) - dr2 = approximate_controls(model, order=2) - dr3 = approximate_controls(model, order=3) - - assert(dr1.order==1) - assert(dr1.X_s.ndim==2) - assert(dr3.X_ss.ndim==3) - -if __name__ == '__main__': - test_fga_higher_order_perturbations() diff --git a/trash/dolo/tests/test_global_2.py b/trash/dolo/tests/test_global_2.py deleted file mode 100644 index 0eb8c668..00000000 --- a/trash/dolo/tests/test_global_2.py +++ /dev/null @@ -1,33 +0,0 @@ -def test_time_iteration_smolyak(): - from dolo import yaml_import, time_iteration - - - filename = 'examples/models/compat/rbc.yaml' - - model = yaml_import(filename) - - import time - - dr = time_iteration(model, pert_order=1, maxit=500, verbose=True) - # dr = time_iteration(model, pert_order=1, maxit=5, smolyak_order=5, verbose=True) - - -def test_time_iteration_spline(): - - import time - from dolo import yaml_import, time_iteration - - - filename = 'examples/models/compat/rbc.yaml' - - model = yaml_import(filename) - print(model.__class__) - - - dr = time_iteration(model, pert_order=1, maxit=5, verbose=True) - # dr = time_iteration(model, pert_order=1, maxit=5, verbose=True) - - -if __name__ == '__main__': - test_time_iteration_spline() - test_time_iteration_smolyak() diff --git a/trash/dolo/tests/test_global_msplines.py b/trash/dolo/tests/test_global_msplines.py deleted file mode 100644 index 1c0e16a7..00000000 --- a/trash/dolo/tests/test_global_msplines.py +++ /dev/null @@ -1,19 +0,0 @@ -import unittest - - -class TestGlobal(unittest.TestCase): - - def test_global_solution(self): - - from dolo import yaml_import, time_iteration - - - filename = 'examples/models/compat/rbc.yaml' - - model = yaml_import(filename) - - import time - - t1 = time.time() - - dr = time_iteration(model, pert_order=1, maxit=5, verbose=True) diff --git a/trash/dolo/tests/test_import_dynare.py b/trash/dolo/tests/test_import_dynare.py deleted file mode 100644 index c4c5affa..00000000 --- a/trash/dolo/tests/test_import_dynare.py +++ /dev/null @@ -1,23 +0,0 @@ - -from dolo.compiler.import_dynare import import_dynare -import os - - -def test_dynare_import(): - - examples_dir = "examples/dynare_modfiles/" - listdir = os.listdir(examples_dir) - - failed = [] - for fname in listdir: - try: - model = import_dynare(examples_dir + fname) - except Exception as e: - failed.append((fname, e)) - - print("Failed:") - for f in failed: print(f) - - -if __name__ == '__main__': - test_dynare_import() diff --git a/trash/dolo/tests/test_mfg_model.py b/trash/dolo/tests/test_mfg_model.py deleted file mode 100644 index 08024474..00000000 --- a/trash/dolo/tests/test_mfg_model.py +++ /dev/null @@ -1,39 +0,0 @@ -def test_mfg_model(): - - from dolo import yaml_import - - # from dolo.algos.commands import time_iteration, simulate, evaluate_policy - from dolo.algos.dtmscc.time_iteration import time_iteration - from dolo.algos.dtmscc.value_iteration import evaluate_policy - from dolo.algos.dtmscc.simulations import simulate, plot_decision_rule - - - model = yaml_import("examples/models/compat/sudden_stop.yaml") - print(model.exogenous) - print(model.is_dtmscc()) - mdr = time_iteration(model) - - # - sim = simulate(model, mdr, 0, horizon=50, n_exp=0) # irf - assert(sim.shape==(50,6)) - sim = simulate(model, mdr, 0, horizon=50, n_exp=1) # one stochastic simulation - assert(sim.shape==(50,6)) - sim = simulate(model, mdr, 0, horizon=50, n_exp=10) # many stochastic simulations - assert(sim.shape==(10,50,6)) - - - mv = evaluate_policy(model, mdr, verbose=True, maxit=500) - - mv = evaluate_policy(model, mdr, verbose=True, maxit=500, initial_guess=mv) - - sim_v = simulate(model, mdr, 0, drv=mv, horizon=50, n_exp=10) - sim_irf = simulate(model, mdr, 0, drv=mv, horizon=50, n_exp=0) - - sim_dr = plot_decision_rule(model, mdr, 'l') - print(sim_dr) - - - -if __name__ == '__main__': - - test_mfg_model() diff --git a/trash/dolo/tests/test_modfile_import.py b/trash/dolo/tests/test_modfile_import.py deleted file mode 100644 index d12d9d9f..00000000 --- a/trash/dolo/tests/test_modfile_import.py +++ /dev/null @@ -1,67 +0,0 @@ -# To change this template, choose Tools | Templates -# and open the template in the editor. - -import unittest -import os - -from trash.dolo.misc.modfile import parse_dynare_text - - -DYNARE_MODFILES_PATH = 'examples/dynare_modfiles/' -exclude = [ - '.directory', # not a modfile - - 't_lag2_check.mod', # 'periods=20;' is not understood - 't_lag2_checka.mod', # 'periods=20;' is not understood - 't_lag2a.mod', # 'periods=20;' is not understood - 't_lag2b.mod', # 'periods=20;' is not understood - - - 'variance_0.mod', # 'periods=20;' is not understood - 'test_matlab.mod', # 'periods=20;' is not understood - - - 'ramsey.mod', # 'old-style shocks' - 'osr_example.mod', # 'old-style shocks' - - 'ramst.mod', # 'deterministic shocks' - 'ramst_a.mod', # 'deterministic shocks' - 'ramst2.mod', # 'deterministic shocks' - - - 't_periods_a.mod', # 'deterministic shocks' - 'example1_varexo_det.mod', # 'deterministic shocks' - 'ramst_normcdf.mod', # uses external function - - 'predetermined_variables.mod', # predetermined variables - -] - -class DynareModfileImportTestCase(unittest.TestCase): - - - def test_dynare_modfile_import(self): - # we test whether each of the modfile included with dynare - # can be imported successfully - - modfiles = os.listdir(DYNARE_MODFILES_PATH) - results = [] - for mf in [m for m in modfiles if m not in exclude]: - fname = DYNARE_MODFILES_PATH + mf - try: - f = file(fname) - txt = f.read() - f.close() - model = parse_dynare_text(txt) - res = (fname,'ok') - results.append( res ) - except Exception as e: - res = (fname, e) - results.append( res ) - print(res) - for r in results: - print(r) - -if __name__ == '__main__': - unittest.main() - diff --git a/trash/dolo/tests/test_parameterized_expectations.py b/trash/dolo/tests/test_parameterized_expectations.py deleted file mode 100644 index 14b481bd..00000000 --- a/trash/dolo/tests/test_parameterized_expectations.py +++ /dev/null @@ -1,38 +0,0 @@ -def test_parameterized_expectations_direct(): - - from dolo import yaml_import - from dolo.algos.dtcscc.parameterized_expectations import parameterized_expectations - from dolo.algos.dtcscc.time_iteration import time_iteration_direct - - model = yaml_import("examples/models/compat/rbc_full.yaml") - - dr_ti = time_iteration_direct(model) - dr_pea = parameterized_expectations(model, direct=True) - - x_ti = dr_ti(dr_ti.grid) - x_pea = dr_pea(dr_ti.grid) - - print(abs(x_ti - x_pea).max() < 1e-5) - -if __name__ == '__main__': - test_parameterized_expectations_direct() - - -def test_parameterized_expectations(): - - from dolo import yaml_import - from dolo.algos.dtcscc.parameterized_expectations import parameterized_expectations - from dolo.algos.dtcscc.time_iteration import time_iteration - - model = yaml_import("examples/models/compat/rbc_full.yaml") - - dr_ti = time_iteration(model) - dr_pea = parameterized_expectations(model, direct=False) - - x_ti = dr_ti(dr_ti.grid) - x_pea = dr_pea(dr_ti.grid) - - print(abs(x_ti - x_pea).max() < 1e-5) - -if __name__ == '__main__': - test_parameterized_expectations() diff --git a/trash/dolo/tests/test_perturbation.py b/trash/dolo/tests/test_perturbation.py deleted file mode 100644 index c02d9b60..00000000 --- a/trash/dolo/tests/test_perturbation.py +++ /dev/null @@ -1,42 +0,0 @@ -import unittest - -class TestPerturbationCase(unittest.TestCase): - - def test_modfile_solve(self): - - from dolo import solve_decision_rule, dynare_import - - model = dynare_import('examples/dynare_modfiles/example1.mod') - - dr = solve_decision_rule(model, order=1) - - # at first order the decision rule is such that: - # y_{t} - ybar = A (y_{t-1} - ybar) + B e_t - - - - def test_dynare_compatibility(self): - - from dolo import solve_decision_rule, dynare_import - - model = dynare_import('examples/dynare_modfiles/example1.mod') - - dr = solve_decision_rule(model, order=1) - - # at first order the decision rule is such that: - # y_{t} - ybar = A (y_{t-1} - ybar) + B e_t - - from dolo.numeric.decision_rules import DynareDecisionRule - dr = DynareDecisionRule(dr, model) - print(dr['ys']) - print(dr['g_a']) - print(dr['g_e']) - - # it can be compared directly to dynare's decision rules - # warning: Dynare's ordering is special - print(dr.ghx) - print(dr.ghu) - -if __name__ == '__main__': - - unittest.main() \ No newline at end of file diff --git a/trash/dolo/tests/test_pickling.py b/trash/dolo/tests/test_pickling.py deleted file mode 100644 index ed37d7a3..00000000 --- a/trash/dolo/tests/test_pickling.py +++ /dev/null @@ -1,32 +0,0 @@ -import unittest -import pickle - -from trash.dolo.symbolic.symbolic import TSymbol - -class PicklingTestCase(unittest.TestCase): - - def test_pickle_variable(self): - - v = TSymbol('v') - eq = v + v(1) + v(2) - - to_be_pickled = { - 'v': v(1), - 'eq': eq - } - - save_string = pickle.dumps(to_be_pickled) - - loaded = pickle.loads(save_string) - - v = loaded['v'] - print(v.__class__) - print(v.date) - print( (v + v)**2 ) - - eqq = loaded['eq'] - print(eqq ) - - -if __name__ == '__main__': - unittest.main() \ No newline at end of file diff --git a/trash/dolo/tests/test_rbc.py b/trash/dolo/tests/test_rbc.py deleted file mode 100644 index 47e9942b..00000000 --- a/trash/dolo/tests/test_rbc.py +++ /dev/null @@ -1,33 +0,0 @@ -from __future__ import print_function - -from dolo import * -import time - -def test_rbc_model(): - - - model = yaml_import('examples/models/compat/rbc.yaml') - - dr = approximate_controls(model) - drg = time_iteration(model) - t1 = time.time() - drg = time_iteration(model) - t2 = time.time() - - sim_dr = plot_decision_rule(model,dr,'k') - - from dolo.algos.dtcscc.value_iteration import evaluate_policy - - pol = evaluate_policy(model, dr, verbose=True) - polg = evaluate_policy(model, drg, verbose=True) - - sim = simulate(model, dr, n_exp=0) # irf - sim = simulate(model, dr, n_exp=2) # stochastic simulations (2 draws) - # extract first simulation - assert(len(sim[0]['k'])==40) - - - -if __name__ == '__main__': - - test_rbc_model() diff --git a/trash/dolo/tests/test_simulations.py b/trash/dolo/tests/test_simulations.py deleted file mode 100644 index e5eb32ef..00000000 --- a/trash/dolo/tests/test_simulations.py +++ /dev/null @@ -1,34 +0,0 @@ -from dolo import * - -def test_simulations_dtcscc(): - model = yaml_import("examples/models/compat/rbc.yaml") - dr = approximate_controls(model) - sim = plot_decision_rule(model,dr,'k') - -def test_simulations_dtcscc_getcorrectbounds(): - model = yaml_import("examples/models/compat/rbc.yaml") - kss = model.get_calibration('k') - - dr = approximate_controls(model) - sim = plot_decision_rule(model,dr,'k',bounds=[1.0,10.0]) - assert(sim['k'].iloc[0] == 1.0) - assert(sim['k'].iloc[-1] == 10.0) - -def test_simulations_dtcscc_getcorrectbounds2(): - model = yaml_import("examples/models/compat/rbc.yaml") - dr = approximate_controls(model) - sim = plot_decision_rule(model,dr,'k') - assert(not(sim['k'].iloc[0] == 1.0)) - assert(not(sim['k'].iloc[-1] == 10.0)) - -def test_simulations_dtmscc(): - model = yaml_import("examples/models/compat/rbc_dtmscc.yaml") - dr = time_iteration(model) - sim = plot_decision_rule(model,dr,'k') - -def test_simulations_dtmscc_getcorrectbounds(): - model = yaml_import("examples/models/compat/rbc_dtmscc.yaml") - dr = time_iteration(model) - sim = plot_decision_rule(model,dr,'k',bounds=[1.0,10.0]) - assert(sim['k'].iloc[0] == 1.0) - assert(sim['k'].iloc[-1] == 10.0) diff --git a/trash/dolo/tests/test_states_perturbation.py b/trash/dolo/tests/test_states_perturbation.py deleted file mode 100644 index 77926e9f..00000000 --- a/trash/dolo/tests/test_states_perturbation.py +++ /dev/null @@ -1,67 +0,0 @@ -import unittest - -class StatesPerturbationsTestCase(unittest.TestCase): - - def test_second_order_accuracy(self): - - # This solves the optimal growth example at second order - # and computes the second order correction to the steady-state - # We test that both the statefree method and the perturbation to states - # yield the same result. - - from dolo.algos.perturbations import yaml_import - model = yaml_import('examples/models/compat/rbc.yaml', compiler=None) - - - from trash.dolo.numeric.perturbations import solve_decision_rule - from trash.dolo.numeric.perturbations_to_states import approximate_controls - - - coeffs = approximate_controls(model,order=2, return_dr=False) - state_perturb = coeffs[0] - - dr = solve_decision_rule(model) - statefree_perturb = dr['ys'] + dr['g_ss']/2.0 - ctls = model.symbols_s['controls'] - ctls_ind = [model.variables.index(v) for v in ctls] - - # the two methods should yield exactly the same result - - from numpy.testing import assert_almost_equal - A = statefree_perturb[ctls_ind] - B = state_perturb - - assert_almost_equal(A, B) # we compare the risk-adjusted constants - - # def test_third_order_accuracy(self): - - - - def test_perturbation_1(self): - from dolo import yaml_import - model = yaml_import('examples/models/compat/rbc.yaml') - from dolo.algos.perturbations import approximate_controls - dr = approximate_controls(model) - - - - def test_perturbation_1_old(self): - from dolo import yaml_import - model = yaml_import('examples/models/compat/rbc.yaml') - from dolo.algos.perturbations import approximate_controls - dr = approximate_controls(model,order=1) - - def test_perturbation_2(self): - from dolo import yaml_import - from dolo.algos.perturbations import approximate_controls - model = yaml_import('examples/models/compat/rbc.yaml') - dr = approximate_controls(model,order=2) - - def test_perturbation_3(self): - from dolo import yaml_import - from dolo.algos.perturbations import approximate_controls - model = yaml_import('examples/models/compat/rbc.yaml') - dr = approximate_controls(model,order=3) - -if __name__ == '__main__': - unittest.main() diff --git a/trash/temp_test_iid_processes.py b/trash/temp_test_iid_processes.py deleted file mode 100644 index 123643a1..00000000 --- a/trash/temp_test_iid_processes.py +++ /dev/null @@ -1,391 +0,0 @@ -from matplotlib import pyplot as plt -import scipy -from scipy.integrate import quad -import scipy.special as special -import numpy as np -from dolo.numeric.processes_iid import * - -σ = 0.1 -μ = 0.2 -N = 10 - -x=2 - -assert(x==2) - -## Polynomial -def f(x): - return x**2 - -## Discretize Normal distribution with GH and EP methods - -distNorm = UNormal(mu=μ, sigma=σ) -disNorm_gh = distNorm.discretize() -disNorm_ep = distNorm.discretize(N=N, method='equiprobable') - -disNorm_gh.n_inodes(0) # number of nodes -# disNorm_gh.inode(0,4) # 4th node -# disNorm_gh.iweight(0,4) # weight associated to 4th node - -######################################################################## -############################# Plots #################################### -######################################################################## -## Plot Equiprobable -weights_ep, nodes_ep = np.array( [*zip(*[*disNorm_ep.iteritems(0)])] ) -plt.plot(nodes_ep, nodes_ep*0, '.') -xl = plt.xlim() -xvec = np.linspace(xl[0], xl[1], 100) -pdf = scipy.stats.norm.pdf(xvec, loc=μ, scale=σ) -plt.plot(xvec, pdf) -plt.grid() - -## Plot Gauss Hermite nodes -weights_gh, nodes_gh = np.array( [*zip(*[*disNorm_gh.iteritems(0)])] ) -plt.plot(nodes_gh, nodes_gh*0, '.') -xl = plt.xlim() -xvec = np.linspace(xl[0], xl[1], 100) -pdf = scipy.stats.norm.pdf(xvec, loc=μ, scale=σ) -plt.plot(xvec, pdf) -plt.grid() - -## Plot Random draws -M=1000 -s_MC = np.random.normal(μ, σ, M) - -count, bins, ignored = plt.hist(s_MC, 30, density=True) -plt.plot(bins, 1/(σ * np.sqrt(2 * np.pi)) * np.exp( - (bins - μ)**2 / (2 * σ **2) ), linewidth=2, color='r') -plt.show() - - -######################################################################## - -## Compute expectation of the normally distributed variable - -expval_normal = quad(lambda x: f(x)/np.sqrt(2*np.pi*σ**2)*np.exp(-(x-μ)**2/(2*σ**2)), -np.inf,np.inf)[0] - -## Compute the mean of random draws -expval_MC = np.array([f(s_MC[j]) for j in range(0,M)]).sum() / M - -## Compute ∑(f(ϵ)*w_ϵ) for each discretization method - -expval_gh = np.array([f(disNorm_gh.inode(0,j))*disNorm_gh.iweight(0,j) for j in range(disNorm_gh.n_inodes(0))]).sum() -expval_ep = np.array([f(disNorm_ep.inode(0,j))*disNorm_ep.iweight(0,j) for j in range(disNorm_ep.n_inodes(0))]).sum() - -## Compare - -expval_normal -expval_gh -expval_ep -expval_MC - -assert(abs(expval_gh-expval_normal)<0.1) -assert(abs(expval_ep-expval_normal)<0.1) - - -################################################################## -############################# UNIFORM ########################### -################################################################## - -a = -1 -b = 1 #### CANNOT PUT 0 -distUni = Uniform(a, b) -disUni = distUni.discretize(N=10) -### Here there is an issue with the uniform - -## Random draws -M=1000 -s_MC = np.random.uniform(a, b, M) - -count, bins, ignored = plt.hist(s_MC, 10, density=True) -plt.plot(bins, np.ones_like(bins)*(1/(b-a)), linewidth=2, color='r') -plt.show() - -## Plot Equiprobable - -weights_uni, nodes_uni = np.array( [*zip(*[*disUni.iteritems(0)])] ) -plt.plot(nodes_uni, nodes_uni*0, '.') -xl = plt.xlim(a,b) -xvec = np.linspace(xl[0], xl[1], 100) -pdf = scipy.stats.uniform.pdf(xvec,a,b-a) -plt.plot(xvec, pdf) -plt.grid() - -## Compute the mean of random draws -expval_MC = np.array([f(s_MC[j]) for j in range(0,M)]).sum() / M -expval_MC - - -# Compute ∑(f(ϵ)*w_ϵ) for each discretization method -expval_ep = np.array([f(disUni.inode(0,j))*disUni.iweight(0,j) for j in range(disUni.n_inodes(0))]).sum() - -expval_ep -expval_MC - -################################################################## -############################# LOGNORMAL ########################### -################################################################## - -σ = 0.1 -μ = 0.3 - -logn = LogNormal(μ=μ, σ=σ) -logn.μ -logn.σ - - -distLog = LogNormal(μ, σ) -disLog = distLog.discretize(N=10) -### Here there is an issue with the uniform - - -## Random draws -M=1000 -s_MC = np.random.lognormal(μ, σ, M) - -count, bins, ignored = plt.hist(s_MC, 10, density=True) -x = np.linspace(min(bins), max(bins), 10000) -pdf = (np.exp(-(np.log(x) - μ)**2 / (2 * σ**2)) / (x * σ * np.sqrt(2 * np.pi))) -plt.plot(x, pdf, linewidth=2, color='r') -plt.axis('tight') -plt.show() - -## Plot Equiprobable - -weights_Log, nodes_Log = np.array( [*zip(*[*disLog.iteritems(0)])] ) -plt.plot(nodes_Log, nodes_Log*0, '.') -xl = plt.xlim() -xvec = np.linspace(xl[0], xl[1], 100) -xvec -pdf = scipy.stats.lognorm.pdf(xvec,s=σ, loc=μ) -plt.plot(xvec, pdf) -pdf2 = (np.exp(-(np.log(x) - μ)**2 / (2 * σ**2)) / (x * σ * np.sqrt(2 * np.pi))) -plt.plot(x, pdf2) -plt.grid() - - -xvec = np.linspace(0,1,100) -ppf0 = scipy.stats.lognorm.ppf(xvec,s=σ, loc=μ, scale=np.exp(μ) ) -ppf1 = distLog.ppf(xvec) - -plt.plot(xvec, ppf0) -plt.plot(xvec, ppf1) - - - - -## Compute the mean of random draws -expval_MC = np.array([f(s_MC[j]) for j in range(0,M)]).sum() / M - -# Compute ∑(f(ϵ)*w_ϵ) for each discretization method -expval_ep = np.array([f(disLog.inode(0,j))*disLog.iweight(0,j) for j in range(disLog.n_inodes(0))]).sum() - -expval_ep -expval_MC - - - - - - - - -################################################################## -############################# BETA ########################### -################################################################## - - - -α = 2 -β = 5 - -distbeta = Beta(α, β) -disbeta = distbeta.discretize(N=10) -### Here there is an issue with the uniform - -## Random draws -M=1000 -s_MC = np.random.beta(α, β, M) - -count, bins, ignored = plt.hist(s_MC, 10, density=True) -x = np.linspace(min(bins), max(bins), 10000) -pdf = scipy.stats.beta.pdf(x,α, β) -plt.plot(x, pdf, linewidth=2, color='r') -plt.axis('tight') -plt.show() - -## Plot Equiprobable -weights_beta, nodes_beta = np.array( [*zip(*[*disbeta.iteritems(0)])] ) -plt.plot(nodes_beta, nodes_beta*0, '.') -xl = plt.xlim() -xvec = np.linspace(xl[0], xl[1], 100) -pdf = scipy.stats.beta.pdf(xvec,α, β) -plt.plot(xvec, pdf) -plt.grid() - - -## Compute the mean of random draws -expval_MC = np.array([f(s_MC[j]) for j in range(0,M)]).sum() / M - -# Compute ∑(f(ϵ)*w_ϵ) for each discretization method -expval_ep = np.array([f(disbeta.inode(0,j))*disbeta.iweight(0,j) for j in range(disbeta.n_inodes(0))]).sum() - -expval_ep -expval_MC - - - - - -################################################################################# -############################ SOME EXTRA STUFF ################################## -################################################################################# -################################################################################# - - - -logn = LogNormal(μ=μ, σ=σ) -logn.μ -logn.σ - -dp = logn.discretize(N=10) -dp2 = logn.discretize(N=10, method='equiprobable') - - -nodes = np.array([x for (w,x) in dp.iteritems(0)]) - -plt.plot(nodes, nodes*0,'.') -plt.xlim(-1,2) -plt.grid() - - - -res_gh = n.discretize(10) -res_ep = n.discretize(10, method='equiprobable') - -for (w,x) in res_ep.iteritems(0): - print(w,x) - - -# neglect integration nodes whose probability is smaller than 1e-5 -for (w,x) in res_gh.iteritems(0,eps=1e-5): - print(w,x) - -def f(x): - return x**2 - - -val = quad(lambda u: f(u)/np.sqrt(2*np.pi*σ**2)*np.exp(-(u-μ)**2/(2*σ**2)), -5, 5) - - -v0 = sum([f(x)*w for (w,x) in res_gh.iteritems(0)]) -v1 = sum([f(x)*w for (w,x) in res_ep.iteritems(0)]) - -print(v0, v1, val) - -sim = norm.simulate(10000,2) - -sim.mean() -sim.std() - - - -dis = n.discretize(N=50, method='equiprobable') - -weights, nodes = np.array( [*zip(*[*dis.iteritems(0)])] ) - - -plt.plot(nodes, nodes*0, '.') -xl = plt.xlim() - - -xvec = np.linspace(xl[0], xl[1], 100) -pdf = scipy.stats.norm.pdf(xvec) -plt.plot(xvec, pdf) -plt.grid() - - -from scipy.stats import beta -import matplotlib.pyplot as plt -import numpy as np -a = 2 -b = 2 -x = np.arange (-50, 50, 0.1) -y = beta.pdf(x,a,b, scale=100, loc=-50) -plt.plot(x,y) - - - -from scipy.stats import beta -import matplotlib.pyplot as plt -import numpy as np -a = 2 -b = 4 -x = np.linspace(-50, 50, 100) -y1 = beta.pdf(x,2, 2, scale=100, loc=-50) -y2 = beta.pdf(x,3, 3, scale=100, loc=-50) -y3 = beta.pdf(x, 4,4, scale=100, loc=-50) - -plt.plot(x, y1, "-", x, y2, "r--", x, y3, "g--") - - - -########################### -logn = LogNormal(μ=μ, σ=σ) -logn.μ -logn.σ - -dp = logn.discretize(N=10) -dp2 = logn.discretize(N=10, method='equiprobable') - - -nodes = np.array([x for (w,x) in dp.iteritems(0)]) - -plt.plot(nodes, nodes*0,'.') -plt.xlim(-1,2) -plt.grid() - - - -res_gh = n.discretize(10) -res_ep = n.discretize(10, method='equiprobable') - -for (w,x) in res_ep.iteritems(0): - print(w,x) - - -# neglect integration nodes whose probability is smaller than 1e-5 -for (w,x) in res_gh.iteritems(0,eps=1e-5): - print(w,x) - -def f(x): - return x**2 - - -val = quad(lambda u: f(u)/np.sqrt(2*np.pi*σ**2)*np.exp(-(u-μ)**2/(2*σ**2)), -5, 5) - - -v0 = sum([f(x)*w for (w,x) in res_gh.iteritems(0)]) -v1 = sum([f(x)*w for (w,x) in res_ep.iteritems(0)]) - -print(v0, v1, val) - -sim = norm.simulate(10000,2) - -sim.mean() -sim.std() - - - -dis = n.discretize(N=50, method='equiprobable') - -weights, nodes = np.array( [*zip(*[*dis.iteritems(0)])] ) - - -plt.plot(nodes, nodes*0, '.') -xl = plt.xlim() - - -xvec = np.linspace(xl[0], xl[1], 100) -pdf = scipy.stats.norm.pdf(xvec) -plt.plot(xvec, pdf) -plt.grid()