From 2b136cb97a257f905d5e1e69c5381f25a933058b Mon Sep 17 00:00:00 2001
From: Daniel Sali <daniel.sali@alphacruncher.com>
Date: Fri, 30 May 2025 16:42:10 +0200
Subject: [PATCH 1/6] Changes following code review
---
src/dynare/__init__.py | 41 ++-
src/dynare/distributions.py | 200 ------------
src/dynare/dynare.py | 323 +++++++++---------
src/dynare/dynare_context.py | 618 +++++++++++++++++++++++++----------
src/dynare/errors.py | 5 +-
src/dynare/indices.py | 123 +++++++
src/dynare/juliapkg.json | 4 +
src/dynare/workspace.py | 134 ++++----
8 files changed, 821 insertions(+), 627 deletions(-)
delete mode 100644 src/dynare/distributions.py
create mode 100644 src/dynare/indices.py
diff --git a/src/dynare/__init__.py b/src/dynare/__init__.py
index 327b6f5..81b1582 100644
--- a/src/dynare/__init__.py
+++ b/src/dynare/__init__.py
@@ -1,20 +1,47 @@
-import subprocess
import logging
+import os # Added import
from juliacall import Main as jl
from .dynare import dynare
-logging.basicConfig(
- level=logging.INFO, format="%(asctime)s - %(name)s - %(levelname)s - %(message)s"
-)
-
logger = logging.getLogger("DynarePython")
+def setup_logging(level_str: str = "INFO", format_string: str = "%(asctime)s - %(name)s - %(levelname)s - %(message)s"):
+ """
+ Set up logging for the DynarePython library.
+
+ Allows client code to explicitly configure logging level and format.
+ This function configures the root logger. If you have a more complex
+ logging setup, you might want to configure the "DynarePython" logger directly.
+ The logging level can also be set using the DYNARE_PYTHON_LOGLEVEL environment variable.
+
+ Args:
+ level_str (str): The desired logging level (e.g., "DEBUG", "INFO", "WARNING").
+ Defaults to "INFO". Is overridden by DYNARE_PYTHON_LOGLEVEL if set.
+ format_string (str): The logging format string.
+ Defaults to "%(asctime)s - %(name)s - %(levelname)s - %(message)s".
+ """
+ env_level_str = os.getenv("DYNARE_PYTHON_LOGLEVEL")
+ if env_level_str:
+ level_str = env_level_str
+
+ level_map = {
+ "DEBUG": logging.DEBUG,
+ "INFO": logging.INFO,
+ "WARNING": logging.WARNING,
+ "ERROR": logging.ERROR,
+ "CRITICAL": logging.CRITICAL,
+ }
+ log_level = level_map.get(level_str.upper(), logging.INFO)
+ logging.basicConfig(level=log_level, format=format_string)
+ logger.info(f"DynarePython logging configured to level: {level_str.upper()}")
+
+
def check_julia_and_dynare():
# Check Julia
julia_path = jl.seval("Sys.BINDIR")
- julia_project = jl.seval("Base.active_project()")
+ _ = jl.seval("Base.active_project()")
logger.info(f"Using julia at: {julia_path}")
# Installed packages
@@ -35,4 +62,4 @@ def _run_once():
_run_once()
-__all__ = ["dynare"]
\ No newline at end of file
+__all__ = ["dynare", "setup_logging"]
\ No newline at end of file
diff --git a/src/dynare/distributions.py b/src/dynare/distributions.py
deleted file mode 100644
index 3c41e07..0000000
--- a/src/dynare/distributions.py
+++ /dev/null
@@ -1,200 +0,0 @@
-""" Python conversion of the base types from the Julia Distributions package. """
-
-from abc import ABC, abstractmethod
-from typing import Type, Any, Tuple, Union, Self
-from math import prod
-
-## VariateForm and its Subtypes
-
-
-class VariateForm(ABC):
- """Specifies the form or shape of the variate or a sample."""
-
- pass
-
-
-class ArrayLikeVariate(VariateForm):
- """Specifies the number of axes of a variate or a sample."""
-
- def __init__(self, n: int):
- self.n = n
-
-
-class Univariate(ArrayLikeVariate):
- def __init__(self):
- super().__init__(0)
-
-
-class Multivariate(ArrayLikeVariate):
- def __init__(self):
- super().__init__(1)
-
-
-class Matrixvariate(ArrayLikeVariate):
- def __init__(self):
- super().__init__(2)
-
-
-class CholeskyVariate(VariateForm):
- """Specifies that the variate or sample is of type Cholesky."""
-
- pass
-
-
-## ValueSupport and its Subtypes
-
-
-class ValueSupport(ABC):
- """Abstract type that specifies the support of elements of samples."""
-
- pass
-
-
-class Discrete(ValueSupport):
- """Represents the support of a discrete random variable (countable)."""
-
- pass
-
-
-class Continuous(ValueSupport):
- """Represents the support of a continuous random variable (uncountable)."""
-
- pass
-
-
-# Promotion rule: Combination of Discrete and Continuous yields Continuous
-def promote_rule(type1: Type[Continuous], type2: Type[Discrete]) -> Type[Continuous]:
- return Continuous
-
-
-## Sampleable
-
-
-class Sampleable(ABC):
- """Sampleable is any type able to produce random values."""
-
- def __init__(
- self, variate_form: Type[VariateForm], value_support: Type[ValueSupport]
- ):
- self.variate_form = variate_form
- self.value_support = value_support
-
- @abstractmethod
- def rand(self) -> Any:
- """Abstract method to generate a random sample."""
- pass
-
-
-def variate_form(sampleable: Type[Sampleable]) -> Type[VariateForm]:
- return sampleable.variate_form
-
-
-def value_support(sampleable: Type[Sampleable]) -> Type[ValueSupport]:
- return sampleable.value_support
-
-
-# Define the size and length of Sampleable objects
-def length(s: Sampleable) -> int:
- if isinstance(s.variate_form, Univariate):
- return 1
- elif isinstance(s.variate_form, Multivariate):
- raise NotImplementedError("Multivariate length is not implemented")
- else:
- return prod(size(s))
-
-
-def size(s: Sampleable) -> Tuple:
- if isinstance(s.variate_form, Univariate):
- return ()
- elif isinstance(s.variate_form, Multivariate):
- return (length(s),)
- else:
- return tuple()
-
-
-def eltype(sampleable: Type[Sampleable]) -> Type:
- if isinstance(sampleable.value_support, Discrete):
- return int
- elif isinstance(sampleable.value_support, Continuous):
- return float
- else:
- raise TypeError("Unknown value support type")
-
-
-## nsamples
-
-
-def nsamples(sampleable: Type[Sampleable], x: Any) -> int:
- if isinstance(sampleable.variate_form, Univariate):
- if isinstance(x, (int, float)):
- return 1
- elif isinstance(x, (list, tuple)):
- return len(x)
- elif isinstance(sampleable.variate_form, Multivariate):
- if isinstance(x, (list, tuple)):
- return 1
- elif isinstance(x, (list, tuple)) and isinstance(x[0], list):
- return len(x[0])
- raise TypeError(f"Unsupported type for nsamples: {type(x)}")
-
-
-## Equality and Hashing for Sampleable objects
-
-
-def sampleable_equal(s1: Sampleable, s2: Sampleable) -> bool:
- if s1.__class__.__name__ != s2.__class__.__name__:
- return False
- for f in vars(s1):
- if getattr(s1, f) != getattr(s2, f):
- return False
- return True
-
-
-def sampleable_hash(s: Sampleable, h: int) -> int:
- hashed = hash(Sampleable)
- hashed = hash(s.__class__.__name__) + h
- for f in vars(s):
- hashed += hash(getattr(s, f))
- return hashed
-
-
-## Distribution as a subtype of Sampleable
-
-
-class Distribution(Sampleable):
- """Distribution is a Sampleable generating random values from a probability distribution."""
-
- def __init__(
- self, variate_form: Type[VariateForm], value_support: Type[ValueSupport]
- ):
- super().__init__(variate_form, value_support)
-
- @classmethod
- def convert_jl_distribution(cls, jl_distribution: Any) -> Self:
- return Distribution()
-
-
-# Defining common distribution types
-class UnivariateDistribution(Distribution):
- def __init__(self, value_support: Type[ValueSupport]):
- super().__init__(Univariate(), value_support)
-
-
-class MultivariateDistribution(Distribution):
- def __init__(self, value_support: Type[ValueSupport]):
- super().__init__(Multivariate(), value_support)
-
-
-class MatrixDistribution(Distribution):
- def __init__(self, value_support: Type[ValueSupport]):
- super().__init__(Matrixvariate(), value_support)
-
-
-class DiscreteDistribution(Distribution):
- def __init__(self, variate_form: Type[VariateForm]):
- super().__init__(variate_form, Discrete)
-
-
-class ContinuousDistribution(Distribution):
- def __init__(self, variate_form: Type[VariateForm]):
- super().__init__(variate_form, Continuous)
diff --git a/src/dynare/dynare.py b/src/dynare/dynare.py
index bb17fab..aa73fb9 100644
--- a/src/dynare/dynare.py
+++ b/src/dynare/dynare.py
@@ -8,189 +8,192 @@ from .errors import DynareError
logger = logging.getLogger("dynare.dynare")
-def dynare(model: str | Path) -> Context:
+def dynare(model: str | Path, *dynare_options: str) -> Context:
try:
if isinstance(model, str):
model = Path(model)
jl.seval("using Serialization")
jl.seval("using Dynare")
- # Double-escape '\' as julia will also interpret the string. Only relevant on Windows
- resolved_path = str(model.resolve()).replace("\\", "\\\\")
+
+ resolved_path_str = str(model.resolve())
+ escaped_path_for_julia_literal = resolved_path_str.replace('\\', '\\\\').replace('"', '\\"')
- jl.seval(
- "using DataFrames, Tables, PythonCall, Dynare, LinearRationalExpectations"
- )
+ julia_options_str_parts = []
+ for opt_str in dynare_options:
+ escaped_opt_for_julia_literal = opt_str.replace('\\', '\\\\').replace('"', '\\"')
+ julia_options_str_parts.append(f'"{escaped_opt_for_julia_literal}"')
+
+ julia_options_cmd_part = " ".join(julia_options_str_parts)
- # Convert the Julia AxisArrayTable fields of a Context to a Pandas DataFrame with PythonCall.pytable
jl.seval(
- """
- using Dynare: EstimationResults, EstimatedParameters, SymbolTable, ModFileInfo
-
- struct PyWork
- analytical_steadystate_variables::Vector{Int}
- data::Py
- datafile::String
- params::Vector{Float64}
- residuals::Vector{Float64}
- dynamic_variables::Vector{Float64}
- exogenous_variables::Vector{Float64}
- observed_variables::Vector{String}
- Sigma_m::Matrix{Float64}
- jacobian::Matrix{Float64}
- qr_jacobian::Matrix{Float64}
- model_has_trend::Vector{Bool}
- histval::Matrix{Union{Float64,Missing}}
- homotopy_setup::Vector{NamedTuple{(:name, :type, :index, :endvalue, :startvalue), Tuple{Symbol, SymbolType, Int64, Float64, Union{Float64, Missing}}}}
- initval_endogenous::Matrix{Union{Float64,Missing}}
- initval_exogenous::Matrix{Union{Float64,Missing}}
- initval_exogenous_deterministic::Matrix{Union{Float64,Missing}}
- endval_endogenous::Matrix{Union{Float64,Missing}}
- endval_exogenous::Matrix{Union{Float64,Missing}}
- endval_exogenous_deterministic::Matrix{Union{Float64,Missing}}
- scenario::Dict{PeriodsSinceEpoch, Dict{PeriodsSinceEpoch, Dict{Symbol, Pair{Float64, Symbol}}}}
- shocks::Vector{Float64}
- perfect_foresight_setup::Dict{String,Any}
- estimated_parameters::EstimatedParameters
- end
-
- function convert_to_pywork(work::Work)::PyWork
- # Convert the AxisArrayTable data to a Pandas DataFrame using pytable
- py_data = pytable(work.data)
-
- return PyWork(
- work.analytical_steadystate_variables,
- py_data,
- work.datafile,
- work.params,
- work.residuals,
- work.dynamic_variables,
- work.exogenous_variables,
- work.observed_variables,
- work.Sigma_m,
- work.jacobian,
- work.qr_jacobian,
- work.model_has_trend,
- work.histval,
- work.homotopy_setup,
- work.initval_endogenous,
- work.initval_exogenous,
- work.initval_exogenous_deterministic,
- work.endval_endogenous,
- work.endval_exogenous,
- work.endval_exogenous_deterministic,
- work.scenario,
- work.shocks,
- work.perfect_foresight_setup,
- work.estimated_parameters
- )
- end
+ """using DataFrames, Tables, PythonCall, Dynare, LinearRationalExpectations
+ using Dynare: EstimationResults, EstimatedParameters, SymbolTable, ModFileInfo, SymbolType, PeriodsSinceEpoch, Work, Simulation, ModelResults, Trends, Context
- struct PySimulation
- firstperiod::PeriodsSinceEpoch
- lastperiod::PeriodsSinceEpoch
- name::String
- statement::String
- data::Py
- end
+ struct PyWork
+ analytical_steadystate_variables::Vector{Int}
+ data::Py
+ datafile::String
+ params::Vector{Float64}
+ residuals::Vector{Float64}
+ dynamic_variables::Vector{Float64}
+ exogenous_variables::Vector{Float64}
+ observed_variables::Vector{String}
+ Sigma_m::Matrix{Float64}
+ jacobian::Matrix{Float64}
+ qr_jacobian::Matrix{Float64}
+ model_has_trend::Vector{Bool}
+ histval::Matrix{Union{Float64,Missing}}
+ homotopy_setup::Vector{NamedTuple{(:name, :type, :index, :endvalue, :startvalue), Tuple{Symbol, SymbolType, Int64, Float64, Union{Float64, Missing}}}}
+ initval_endogenous::Matrix{Union{Float64,Missing}}
+ initval_exogenous::Matrix{Union{Float64,Missing}}
+ initval_exogenous_deterministic::Matrix{Union{Float64,Missing}}
+ endval_endogenous::Matrix{Union{Float64,Missing}}
+ endval_exogenous::Matrix{Union{Float64,Missing}}
+ endval_exogenous_deterministic::Matrix{Union{Float64,Missing}}
+ scenario::Dict{PeriodsSinceEpoch, Dict{PeriodsSinceEpoch, Dict{Symbol, Pair{Float64, Symbol}}}}
+ shocks::Vector{Float64}
+ perfect_foresight_setup::Dict{String,Any}
+ estimated_parameters::EstimatedParameters
+ end
+
+ function convert_to_pywork(work::Work)::PyWork
+ # Convert the AxisArrayTable data to a Pandas DataFrame using pytable
+ py_data = pytable(work.data)
- function convert_to_pysimulation(simulation::Simulation)::PySimulation
- # Convert the AxisArrayTable data to a Pandas DataFrame using pytable
- py_data = pytable(simulation.data)
-
- return PySimulation(
- simulation.firstperiod,
- simulation.lastperiod,
- simulation.name,
- simulation.statement,
- py_data
- )
- end
-
+ return PyWork(
+ work.analytical_steadystate_variables,
+ py_data,
+ work.datafile,
+ work.params,
+ work.residuals,
+ work.dynamic_variables,
+ work.exogenous_variables,
+ work.observed_variables,
+ work.Sigma_m,
+ work.jacobian,
+ work.qr_jacobian,
+ work.model_has_trend,
+ work.histval,
+ work.homotopy_setup,
+ work.initval_endogenous,
+ work.initval_exogenous,
+ work.initval_exogenous_deterministic,
+ work.endval_endogenous,
+ work.endval_exogenous,
+ work.endval_exogenous_deterministic,
+ work.scenario,
+ work.shocks,
+ work.perfect_foresight_setup,
+ work.estimated_parameters
+ )
+ end
+
+ struct PySimulation
+ firstperiod::PeriodsSinceEpoch
+ lastperiod::PeriodsSinceEpoch
+ name::String
+ statement::String
+ data::Py
+ end
+
+ function convert_to_pysimulation(simulation::Simulation)::PySimulation
+ # Convert the AxisArrayTable data to a Pandas DataFrame using pytable
+ py_data = pytable(simulation.data)
- # Define the PyModelResult structure with Pandas DataFrame fields
- mutable struct PyModelResult
- irfs::Dict{Symbol, Py}
- trends::Trends
- stationary_variables::Vector{Bool}
- estimation::EstimationResults
- filter::Py # Pandas DataFrame
- forecast::Vector{Py} # Vector of Pandas DataFrames
- initial_smoother::Py # Pandas DataFrame
- linearrationalexpectations::LinearRationalExpectationsResults
- simulations::Vector{PySimulation}
- smoother::Py # Pandas DataFrame
- solution_derivatives::Vector{Matrix{Float64}}
- # sparsegrids::SparsegridsResults
- end
+ return PySimulation(
+ simulation.firstperiod,
+ simulation.lastperiod,
+ simulation.name,
+ simulation.statement,
+ py_data
+ )
+ end
- function convert_to_pymodelresult(model_result::ModelResults)::PyModelResult
- py_irfs = Dict{Symbol, Py}()
- for (key, axis_array_table) in model_result.irfs
- py_irfs[key] = pytable(axis_array_table)
+
+ # Define the PyModelResult structure with Pandas DataFrame fields
+ mutable struct PyModelResult
+ irfs::Dict{Symbol, Py}
+ trends::Trends
+ stationary_variables::Vector{Bool}
+ estimation::EstimationResults
+ filter::Py # Pandas DataFrame
+ forecast::Vector{Py} # Vector of Pandas DataFrames
+ initial_smoother::Py # Pandas DataFrame
+ linearrationalexpectations::LinearRationalExpectationsResults
+ simulations::Vector{PySimulation}
+ smoother::Py # Pandas DataFrame
+ solution_derivatives::Vector{Matrix{Float64}}
+ # sparsegrids::SparsegridsResults
end
- py_forecast = [pytable(forecast) for forecast in model_result.forecast]
-
- return PyModelResult(
- py_irfs,
- model_result.trends,
- model_result.stationary_variables,
- model_result.estimation,
- pytable(model_result.filter),
- py_forecast,
- pytable(model_result.initial_smoother),
- model_result.linearrationalexpectations,
- [convert_to_pysimulation(simulation) for simulation in model_result.simulations],
- pytable(model_result.smoother),
- model_result.solution_derivatives,
- # model_result.sparsegrids
- )
- end
+ function convert_to_pymodelresult(model_result::ModelResults)::PyModelResult
+ py_irfs = Dict{Symbol, Py}()
+ for (key, axis_array_table) in model_result.irfs
+ py_irfs[key] = pytable(axis_array_table)
+ end
- struct PyResults
- model_results::Vector{PyModelResult}
- end
+ py_forecast = [pytable(forecast) for forecast in model_result.forecast]
- struct PyContext
- symboltable::SymbolTable
- models::Vector{Model}
- modfileinfo::ModFileInfo
- results::PyResults # Now holds PyModelResult instead of ModelResult
- work::PyWork
- workspaces::Dict
- end
+ return PyModelResult(
+ py_irfs,
+ model_result.trends,
+ model_result.stationary_variables,
+ model_result.estimation,
+ pytable(model_result.filter),
+ py_forecast,
+ pytable(model_result.initial_smoother),
+ model_result.linearrationalexpectations,
+ [convert_to_pysimulation(simulation) for simulation in model_result.simulations],
+ pytable(model_result.smoother),
+ model_result.solution_derivatives,
+ # model_result.sparsegrids
+ )
+ end
- function convert_to_pycontext(ctx::Context)::PyContext
- # Convert each ModelResults in the Context to PyModelResult
- py_model_results = [convert_to_pymodelresult(model_result) for model_result in ctx.results.model_results]
+ struct PyResults
+ model_results::Vector{PyModelResult}
+ end
+
+ struct PyContext
+ symboltable::SymbolTable
+ models::Vector{Model}
+ modfileinfo::ModFileInfo
+ results::PyResults # Now holds PyModelResult instead of ModelResult
+ work::PyWork
+ workspaces::Dict
+ end
- # Create a PyResults structure with the converted PyModelResults
- py_results = PyResults(py_model_results)
+ function convert_to_pycontext(ctx::Context)::PyContext
+ # Convert each ModelResults in the Context to PyModelResult
+ py_model_results = [convert_to_pymodelresult(model_result) for model_result in ctx.results.model_results]
- # Convert the Work structure
- py_work = convert_to_pywork(ctx.work)
+ # Create a PyResults structure with the converted PyModelResults
+ py_results = PyResults(py_model_results)
- # Return a new PyContext with PyResults and PyWork
- return PyContext(
- ctx.symboltable,
- ctx.models,
- ctx.modfileinfo,
- py_results, # PyResults instead of Results
- py_work,
- ctx.workspaces
- )
- end
- """
- )
- context = jl.seval(
- f"""ctx = @dynare "{resolved_path}";
- if !(ctx isa Dynare.Context)
- throw(error("Failed to produce a Dynare context."))
- else
- convert_to_pycontext(ctx)
+ # Convert the Work structure
+ py_work = convert_to_pywork(ctx.work)
+
+ # Return a new PyContext with PyResults and PyWork
+ return PyContext(
+ ctx.symboltable,
+ ctx.models,
+ ctx.modfileinfo,
+ py_results, # PyResults instead of Results
+ py_work,
+ ctx.workspaces
+ )
end
"""
)
+
+ julia_command = f'''ctx = @dynare "{escaped_path_for_julia_literal}" {julia_options_cmd_part};
+ if !(ctx isa Dynare.Context)
+ throw(error("Failed to produce a Dynare context."))
+ else
+ convert_to_pycontext(ctx)
+ end
+ '''
+ context = jl.seval(julia_command)
return Context.from_julia(context)
except JuliaError as e:
raise DynareError.from_julia_error(e)
diff --git a/src/dynare/dynare_context.py b/src/dynare/dynare_context.py
index a6d5d7d..ea24492 100644
--- a/src/dynare/dynare_context.py
+++ b/src/dynare/dynare_context.py
@@ -3,10 +3,12 @@ import ctypes
import os
from enum import Enum
from dataclasses import dataclass, field
-from typing import List, Any, Dict, Tuple, Union, Self, TypeVar, Generic
+from typing import List, Any, Dict, Tuple, Union
import numpy as np
import pandas as pd
import scipy.stats as stats
+from juliacall import Main as jl # Added import
+from .indices import PyIndices
logger = logging.getLogger("dynare.dynare_context")
@@ -20,17 +22,9 @@ class SymbolType(Enum):
DynareFunction = 4
@classmethod
- def from_julia(cls, jl_symboltype) -> Self:
- return SymbolType(
- int(
- str(repr(jl_symboltype))
- .replace("Julia: Endogenous::SymbolType = ", "")
- .replace("Julia: Exogenous::SymbolType = ", "")
- .replace("Julia: ExogenousDeterministic::SymbolType = ", "")
- .replace("Julia: Parameter::SymbolType = ", "")
- .replace("Julia: DynareFunction::SymbolType = ", "")
- )
- )
+ def from_julia(cls, jl_symboltype: Any) -> "SymbolType":
+ """Converts a Julia SymbolType enum to its Python equivalent."""
+ return cls[str(jl_symboltype)]
class EstimatedParameterType(Enum):
@@ -43,19 +37,9 @@ class EstimatedParameterType(Enum):
EstCorrMeasurement = 6
@classmethod
- def from_julia(cls, jl_estimatedparametertype) -> Self:
- return EstimatedParameterType(
- int(
- str(repr(jl_estimatedparametertype))
- .replace("Julia: EstParameter::EstimatedParameterType = ", "")
- .replace("Julia: EstSDShock::EstimatedParameterType = ", "")
- .replace("Julia: EstSDMeasurement::EstimatedParameterType = ", "")
- .replace("Julia: EstVarShock::EstimatedParameterType = ", "")
- .replace("Julia: EstVarMeasurement::EstimatedParameterType = ", "")
- .replace("Julia: EstCorrShock::EstimatedParameterType = ", "")
- .replace("Julia: EstCorrMeasurement::EstimatedParameterType = ", "")
- )
- )
+ def from_julia(cls, jl_estimatedparametertype: Any) -> "EstimatedParameterType":
+ """Converts a Julia EstimatedParameterType enum to its Python equivalent."""
+ return cls[str(jl_estimatedparametertype)]
def get_member(obj, member, default=None):
@@ -119,7 +103,8 @@ class ModFileInfo:
)
@classmethod
- def from_julia(cls, jl_modfileinfo) -> Self:
+ def from_julia(cls, jl_modfileinfo: Any) -> "ModFileInfo":
+ """Converts a Julia ModFileInfo object to its Python equivalent."""
return ModFileInfo(
endval_is_reset=jl_modfileinfo.endval_is_reset,
has_auxiliary_variables=jl_modfileinfo.has_auxiliary_variables,
@@ -146,86 +131,6 @@ class ModFileInfo:
)
-@dataclass
-class PyIndices:
- current: List[int]
- forward: List[int]
- purely_forward: List[int]
- backward: List[int]
- both: List[int]
- non_backward: List[int]
-
- static: List[int]
- dynamic: List[int]
- dynamic_current: List[int]
-
- current_in_dynamic: List[int]
- forward_in_dynamic: List[int]
- backward_in_dynamic: List[int]
-
- current_in_dynamic_jacobian: List[int]
- current_in_static_jacobian: List[int]
-
- exogenous: List[int]
- n_endogenous: int
-
- D_columns: Tuple[List[int], List[int]]
- E_columns: Tuple[List[int], List[int]]
- UD_columns: List[int]
- UE_columns: List[int]
-
- def __repr__(self) -> str:
- return (
- f"PyIndices(\n"
- f" current={self.current},\n"
- f" forward={self.forward},\n"
- f" purely_forward={self.purely_forward},\n"
- f" backward={self.backward},\n"
- f" both={self.both},\n"
- f" non_backward={self.non_backward},\n"
- f" static={self.static},\n"
- f" dynamic={self.dynamic},\n"
- f" dynamic_current={self.dynamic_current},\n"
- f" current_in_dynamic={self.current_in_dynamic},\n"
- f" forward_in_dynamic={self.forward_in_dynamic},\n"
- f" backward_in_dynamic={self.backward_in_dynamic},\n"
- f" current_in_dynamic_jacobian={self.current_in_dynamic_jacobian},\n"
- f" current_in_static_jacobian={self.current_in_static_jacobian},\n"
- f" exogenous={self.exogenous},\n"
- f" n_endogenous={self.n_endogenous},\n"
- f" D_columns={self.D_columns},\n"
- f" E_columns={self.E_columns},\n"
- f" UD_columns={self.UD_columns},\n"
- f" UE_columns={self.UE_columns}\n"
- f")"
- )
-
- @classmethod
- def from_julia(cls, julia_obj):
- return cls(
- current=list(julia_obj.current),
- forward=list(julia_obj.forward),
- purely_forward=list(julia_obj.purely_forward),
- backward=list(julia_obj.backward),
- both=list(julia_obj.both),
- non_backward=list(julia_obj.non_backward),
- static=list(julia_obj.static),
- dynamic=list(julia_obj.dynamic),
- dynamic_current=list(julia_obj.dynamic_current),
- current_in_dynamic=list(julia_obj.current_in_dynamic),
- forward_in_dynamic=list(julia_obj.forward_in_dynamic),
- backward_in_dynamic=list(julia_obj.backward_in_dynamic),
- current_in_dynamic_jacobian=list(julia_obj.current_in_dynamic_jacobian),
- current_in_static_jacobian=list(julia_obj.current_in_static_jacobian),
- exogenous=list(julia_obj.exogenous),
- n_endogenous=julia_obj.n_endogenous,
- D_columns=(list(julia_obj.D_columns.D), list(julia_obj.D_columns.jacobian)),
- E_columns=(list(julia_obj.E_columns.E), list(julia_obj.E_columns.jacobian)),
- UD_columns=list(julia_obj.UD_columns),
- UE_columns=list(julia_obj.UE_columns),
- )
-
-
@dataclass
class Model:
endogenous_nbr: int
@@ -422,7 +327,8 @@ class Model:
)
@classmethod
- def from_julia(cls, jl_model) -> Self:
+ def from_julia(cls, jl_model: Any) -> "Model":
+ """Converts a Julia Model object to its Python equivalent."""
return Model(
endogenous_nbr=get_member(jl_model, "endogenous_nbr"),
exogenous_nbr=get_member(jl_model, "exogenous_nbr"),
@@ -567,7 +473,8 @@ class Simulation:
)
@classmethod
- def from_julia(cls, jl_simulation) -> Self:
+ def from_julia(cls, jl_simulation: Any) -> "Simulation":
+ """Converts a Julia Simulation object to its Python equivalent."""
return Simulation(
firstperiod=jl_simulation.firstperiod,
lastperiod=jl_simulation.lastperiod,
@@ -611,7 +518,8 @@ class Trends:
)
@classmethod
- def from_julia(cls, jl_trends) -> Self:
+ def from_julia(cls, jl_trends: Any) -> "Trends":
+ """Converts a Julia Trends object to its Python equivalent."""
return Trends(
endogenous_steady_state=list(jl_trends.endogenous_steady_state),
endogenous_terminal_steady_state=list(
@@ -662,7 +570,8 @@ class EstimationResults:
)
@classmethod
- def from_julia(cls, jl_estimationresults) -> Self:
+ def from_julia(cls, jl_estimationresults: Any) -> "EstimationResults":
+ """Converts a Julia EstimationResults object to its Python equivalent."""
return EstimationResults(
posterior_mode=list(jl_estimationresults.posterior_mode),
transformed_posterior_mode=list(
@@ -714,7 +623,8 @@ class LinearRationalExpectationsResults:
)
@classmethod
- def from_julia(cls, jl_linearrationalexpectationsresults) -> Self:
+ def from_julia(cls, jl_linearrationalexpectationsresults: Any) -> "LinearRationalExpectationsResults":
+ """Converts a Julia LinearRationalExpectationsResults object to its Python equivalent."""
return LinearRationalExpectationsResults(
eigenvalues=jl_linearrationalexpectationsresults.eigenvalues.to_numpy(),
g1=jl_linearrationalexpectationsresults.g1.to_numpy(),
@@ -735,15 +645,9 @@ class PySGSolver(Enum):
PATHSolver = 3
@classmethod
- def from_julia(cls, jl_sgsolver) -> Self:
- return PySGSolver(
- int(
- str(repr(jl_sgsolver))
- .replace("Julia: NLsolver::SGSolver = ", "")
- .replace("Julia: NonlinearSolver::SGSolver = ", "")
- .replace("Julia: PATHSolver::SGSolver = ", "")
- )
- )
+ def from_julia(cls, jl_sgsolver: Any) -> "PySGSolver":
+ """Converts a Julia SGSolver enum to its Python equivalent PySGSolver."""
+ return cls[str(jl_sgsolver)]
libtas_path = os.getenv("LIBTASMANIAN_PATH")
@@ -751,6 +655,9 @@ if libtas_path:
TASlib = ctypes.CDLL(libtas_path) # Load the shared library
else:
TASlib = None
+ logger.warning(
+ "LIBTASMANIAN_PATH environment variable not set. Tasmanian-dependent features may not be available."
+ )
@dataclass
@@ -773,7 +680,8 @@ class PyTasmanianSG:
)
@classmethod
- def from_julia(cls, jl_tasmaniansg) -> Self:
+ def from_julia(cls, jl_tasmaniansg: Any) -> "PyTasmanianSG":
+ """Converts a Julia TasmanianSG object to its Python equivalent PyTasmanianSG."""
if not TASlib:
raise ValueError("It looks like Tasmanian is not installed. Please install from https://github.com/ORNL/TASMANIAN and set LIBTASMANIAN_PATH")
c_func = getattr(TASlib, "tsgConstructTasmanianSparseGrid")
@@ -790,13 +698,31 @@ class PyTasmanianSG:
depth=jl_tasmaniansg.depth,
)
-
-T = TypeVar("T") # Placeholder for concrete_jac
-N = TypeVar("N") # Placeholder for name
+ def __del__(self):
+ if hasattr(self, 'pGrid') and self.pGrid and TASlib:
+ try:
+ # Attempt to get the destructor function from the loaded TASlib
+ destructor_func = getattr(TASlib, "tsgDestructTasmanianSparseGrid")
+ destructor_func.argtypes = [ctypes.POINTER(ctypes.c_void_p)]
+ destructor_func.restype = None # Typically void for destructors
+
+ destructor_func(self.pGrid)
+ logger.debug(f"Tasmanian sparse grid {self.pGrid} destructed successfully.")
+ except AttributeError:
+ # This might happen if TASlib is loaded but doesn't have the expected destructor
+ logger.error(
+ "Tasmanian library is loaded but 'tsgDestructTasmanianSparseGrid' function not found. Cannot free pGrid."
+ )
+ except Exception as e:
+ logger.error(f"Error during destruction of Tasmanian sparse grid {self.pGrid}: {e}")
+ elif hasattr(self, 'pGrid') and self.pGrid and not TASlib:
+ logger.warning(
+ f"Tasmanian sparse grid {self.pGrid} might not be freed because TASlib is not loaded (LIBTASMANIAN_PATH not set or invalid)."
+ )
@dataclass
-class PyGeneralizedFirstOrderAlgorithm(Generic[T, N]):
+class PyGeneralizedFirstOrderAlgorithm:
linesearch: Any
trustregion: Any
descent: Any
@@ -821,8 +747,9 @@ class PyGeneralizedFirstOrderAlgorithm(Generic[T, N]):
@classmethod
def from_julia(
cls,
- jl_generalizedfirstorderalgorithm,
- ):
+ jl_generalizedfirstorderalgorithm: Any,
+ ) -> "PyGeneralizedFirstOrderAlgorithm":
+ """Converts a Julia GeneralizedFirstOrderAlgorithm object to its Python equivalent."""
return cls(
linesearch=jl_generalizedfirstorderalgorithm.linesearch,
trustregion=jl_generalizedfirstorderalgorithm.trustregion,
@@ -880,7 +807,8 @@ class PySparsegridsResults:
)
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PySparsegridsResults":
+ """Converts a Julia SparsegridsResults object to its Python equivalent PySparsegridsResults."""
return cls(
average_error=julia_obj.average_error,
average_iteration_time=julia_obj.average_iteration_time,
@@ -943,7 +871,8 @@ class ModelResults:
)
@classmethod
- def from_julia(cls, jl_modelresults) -> Self:
+ def from_julia(cls, jl_modelresults: Any) -> "ModelResults":
+ """Converts a Julia ModelResults object to its Python equivalent."""
return ModelResults(
irfs={
repr(k).replace("Julia: ", ""): v
@@ -977,10 +906,11 @@ class Results:
model_results: List[ModelResults]
def __repr__(self) -> str:
- return f"Results(\n" f" model_results={self.model_results}\n" ")"
+ return f"Results(\\\\n model_results={self.model_results}\\\\n)"
@classmethod
- def from_julia(cls, jl_results) -> Self:
+ def from_julia(cls, jl_results: Any) -> "Results":
+ """Converts a Julia Results object to its Python equivalent."""
return Results(
model_results=[
ModelResults.from_julia(jlmr) for jlmr in jl_results.model_results
@@ -989,49 +919,380 @@ class Results:
def distribution_from_julia(
- jl_distribution,
+ jl_distribution: Any,
) -> Union[
stats.rv_continuous, stats.rv_discrete, stats._multivariate.multi_rv_generic, None
]:
- conversion_map = {
- "Normal": lambda dist: stats.norm(loc=dist.μ, scale=dist.σ),
- "Exponential": lambda dist: stats.expon(scale=1 / dist.θ),
- "Uniform": lambda dist: stats.uniform(loc=dist.a, scale=dist.b - dist.a),
- "Gamma": lambda dist: stats.gamma(a=dist.α, scale=1 / dist.θ),
- "Beta": lambda dist: stats.beta(a=dist.α, b=dist.β),
- "Binomial": lambda dist: stats.binom(n=dist.n, p=dist.p),
- "Poisson": lambda dist: stats.poisson(mu=dist.λ),
- "Geometric": lambda dist: stats.geom(p=dist.p),
- "NegativeBinomial": lambda dist: stats.nbinom(n=dist.r, p=dist.p),
- "ChiSquared": lambda dist: stats.chi2(df=dist.ν),
- "TDist": lambda dist: stats.t(df=dist.ν),
- "Laplace": lambda dist: stats.laplace(loc=dist.μ, scale=dist.β),
- "Cauchy": lambda dist: stats.cauchy(loc=dist.x0, scale=dist.γ),
- "LogNormal": lambda dist: stats.lognorm(s=dist.σ, scale=np.exp(dist.μ)),
- "Weibull": lambda dist: stats.weibull_min(c=dist.k, scale=dist.λ),
- "Gumbel": lambda dist: stats.gumbel_r(loc=dist.μ, scale=dist.β),
- "MvNormal": lambda dist: stats.multivariate_normal(
- mean=np.array(dist.μ), cov=np.array(dist.Σ)
- ),
- "MvTDist": lambda dist: stats.multivariate_t(
- loc=np.array(dist.μ), shape=np.array(dist.Σ), df=dist.ν
- ),
- "Dirichlet": lambda dist: stats.dirichlet(alpha=np.array(dist.α)),
- "Bernoulli": lambda dist: stats.bernoulli(p=dist.p),
- "Categorical": lambda dist: stats.multinomial(n=1, p=np.array(dist.p)),
- "Multinomial": lambda dist: stats.multinomial(n=dist.n, p=np.array(dist.p)),
- }
-
- # Get the Julia distribution type name
- julia_dist_type = str(type(jl_distribution)).split(".")[-1].strip("'>")
-
- # Convert the Julia distribution to a scipy.stats distribution
- if julia_dist_type in conversion_map:
- return conversion_map[julia_dist_type](jl_distribution)
- else:
- raise ValueError(
- f"No conversion available for Julia distribution type: {julia_dist_type}"
- )
+ """
+ Converts a Julia Distributions.jl distribution object to its equivalent scipy.stats object.
+ Enhanced for robustness, safer Julia interaction, and detailed validation.
+ """
+ if jl_distribution is None:
+ return None
+
+ dist_name_for_log = "UnknownDistribution"
+ try:
+ julia_type_obj = jl.typeof(jl_distribution)
+ dist_name = jl.nameof(julia_type_obj).__str__()
+ dist_name_for_log = dist_name # Update for logging once name is known
+
+ # --- Helper functions for parameter extraction ---
+ def _get_scalar_param(param_name_str: str) -> Union[float, None]:
+ try:
+ jl_val = jl.getproperty(jl_distribution, jl.Symbol(param_name_str))
+ return jl.pyconvert(float, jl_val)
+ except Exception:
+ logger.debug(f"Scalar parameter '{param_name_str}' not found or not convertible for {dist_name}.")
+ return None
+
+ def _get_int_param(param_name_str: str) -> Union[int, None]:
+ try:
+ jl_val = jl.getproperty(jl_distribution, jl.Symbol(param_name_str))
+ return jl.pyconvert(int, jl_val)
+ except Exception:
+ logger.debug(f"Integer parameter '{param_name_str}' not found or not convertible for {dist_name}.")
+ return None
+
+ def _get_numpy_array_param(param_name_str: str) -> Union[np.ndarray, None]:
+ try:
+ jl_val = jl.getproperty(jl_distribution, jl.Symbol(param_name_str))
+ return jl.pyconvert(np.ndarray, jl_val)
+ except Exception:
+ logger.debug(f"NumPy array parameter '{param_name_str}' not found or not convertible for {dist_name}.")
+ return None
+
+ # --- Individual conversion functions for each distribution type ---
+ def _convert_normal():
+ mu = _get_scalar_param("μ")
+ sigma = _get_scalar_param("σ")
+ if mu is not None and sigma is not None:
+ if sigma <= 0:
+ logger.error(f"Normal 'σ' must be positive for {dist_name}, got {sigma}.")
+ return None
+ return stats.norm(loc=mu, scale=sigma)
+ return None
+
+ def _convert_exponential():
+ theta = _get_scalar_param("θ") # rate in Distributions.jl
+ if theta is not None:
+ if theta <= 0:
+ logger.error(f"Exponential 'θ' (rate) must be positive for {dist_name}, got {theta}.")
+ return None
+ return stats.expon(scale=1/theta) # scipy.expon uses scale = 1/rate
+ return None
+
+ def _convert_uniform():
+ a = _get_scalar_param("a")
+ b = _get_scalar_param("b")
+ if a is not None and b is not None:
+ if b < a:
+ logger.error(f"Uniform 'b' must be >= 'a' for {dist_name}, got a={a}, b={b}.")
+ return None
+ return stats.uniform(loc=a, scale=b - a)
+ return None
+
+ def _convert_gamma(): # Distributions.jl: Gamma(shape, scale) = Gamma(α, θ)
+ alpha = _get_scalar_param("α") # shape
+ theta = _get_scalar_param("θ") # scale
+ if alpha is not None and theta is not None:
+ if alpha <= 0 or theta <= 0:
+ logger.error(f"Gamma 'α' (shape) and 'θ' (scale) must be positive for {dist_name}, got α={alpha}, θ={theta}.")
+ return None
+ return stats.gamma(a=alpha, scale=theta) # scipy.stats.gamma(a=shape, scale=scale)
+ return None
+
+ def _convert_beta():
+ alpha = _get_scalar_param("α")
+ beta_p = _get_scalar_param("β")
+ if alpha is not None and beta_p is not None:
+ if alpha <= 0 or beta_p <= 0:
+ logger.error(f"Beta 'α' and 'β' must be positive for {dist_name}, got α={alpha}, β={beta_p}.")
+ return None
+ return stats.beta(a=alpha, b=beta_p)
+ return None
+
+ def _convert_binomial():
+ n = _get_int_param("n")
+ p = _get_scalar_param("p")
+ if n is not None and p is not None:
+ if n < 0:
+ logger.error(f"Binomial 'n' must be non-negative for {dist_name}, got n={n}.")
+ return None
+ if not (0 <= p <= 1):
+ logger.error(f"Binomial 'p' must be in [0,1] for {dist_name}, got p={p}.")
+ return None
+ return stats.binom(n=n, p=p)
+ return None
+
+ def _convert_poisson():
+ lambda_ = _get_scalar_param("λ")
+ if lambda_ is not None:
+ if lambda_ < 0: # scipy allows lambda_ == 0
+ logger.error(f"Poisson 'λ' must be non-negative for {dist_name}, got {lambda_}.")
+ return None
+ return stats.poisson(mu=lambda_)
+ return None
+
+ def _convert_geometric():
+ p = _get_scalar_param("p")
+ if p is not None:
+ if not (0 < p <= 1):
+ logger.error(f"Geometric 'p' must be in (0,1] for {dist_name}, got {p}.")
+ return None
+ return stats.geom(p=p)
+ return None
+
+ def _convert_negative_binomial():
+ r = _get_int_param("r") # number of successes
+ p = _get_scalar_param("p") # probability of success
+ if r is not None and p is not None:
+ if r <= 0:
+ logger.error(f"NegativeBinomial 'r' (successes) must be positive for {dist_name}, got {r}.")
+ return None
+ if not (0 < p <= 1): # p can be 1
+ logger.error(f"NegativeBinomial 'p' must be in (0,1] for {dist_name}, got {p}.")
+ return None
+ return stats.nbinom(n=r, p=p)
+ return None
+
+ def _convert_chi_squared(): # Julia: ChiSquared, Scipy: chi2
+ df = _get_scalar_param("ν") # degrees of freedom
+ if df is not None:
+ if df <= 0:
+ logger.error(f"ChiSquared 'ν' (df) must be positive for {dist_name}, got {df}.")
+ return None
+ return stats.chi2(df=df)
+ return None
+
+ def _convert_tdist():
+ df = _get_scalar_param("ν") # degrees of freedom
+ if df is not None:
+ if df <= 0:
+ logger.error(f"TDist 'ν' (df) must be positive for {dist_name}, got {df}.")
+ return None
+ return stats.t(df=df)
+ return None
+
+ def _convert_laplace():
+ mu = _get_scalar_param("μ") # location
+ beta_scale = _get_scalar_param("β") # scale
+ if mu is not None and beta_scale is not None:
+ if beta_scale <= 0:
+ logger.error(f"Laplace 'β' (scale) must be positive for {dist_name}, got {beta_scale}.")
+ return None
+ return stats.laplace(loc=mu, scale=beta_scale)
+ return None
+
+ def _convert_cauchy():
+ x0 = _get_scalar_param("x0") # location
+ gamma_scale = _get_scalar_param("γ") # scale
+ if x0 is not None and gamma_scale is not None:
+ if gamma_scale <= 0:
+ logger.error(f"Cauchy 'γ' (scale) must be positive for {dist_name}, got {gamma_scale}.")
+ return None
+ return stats.cauchy(loc=x0, scale=gamma_scale)
+ return None
+
+ def _convert_log_normal():
+ mu_log = _get_scalar_param("μ") # log-mean
+ sigma_log = _get_scalar_param("σ") # log-std
+ if mu_log is not None and sigma_log is not None:
+ if sigma_log <= 0:
+ logger.error(f"LogNormal 'σ' (log-std) must be positive for {dist_name}, got {sigma_log}.")
+ return None
+ return stats.lognorm(s=sigma_log, scale=np.exp(mu_log)) # s is shape (log-std), scale is exp(log-mean)
+ return None
+
+ def _convert_weibull():
+ k_shape = _get_scalar_param("k") # shape
+ lambda_scale = _get_scalar_param("λ") # scale
+ if k_shape is not None and lambda_scale is not None:
+ if k_shape <= 0 or lambda_scale <= 0:
+ logger.error(f"Weibull 'k' (shape) and 'λ' (scale) must be positive for {dist_name}, got k={k_shape}, λ={lambda_scale}.")
+ return None
+ return stats.weibull_min(c=k_shape, scale=lambda_scale)
+ return None
+
+ def _convert_gumbel():
+ mu_loc = _get_scalar_param("μ") # location
+ beta_scale = _get_scalar_param("β") # scale
+ if mu_loc is not None and beta_scale is not None:
+ if beta_scale <= 0:
+ logger.error(f"Gumbel 'β' (scale) must be positive for {dist_name}, got {beta_scale}.")
+ return None
+ return stats.gumbel_r(loc=mu_loc, scale=beta_scale)
+ return None
+
+ def _convert_mv_normal():
+ mean_val = _get_numpy_array_param("μ")
+ cov_val = _get_numpy_array_param("Σ")
+ if mean_val is not None and cov_val is not None:
+ if mean_val.ndim != 1:
+ logger.error(f"MvNormal 'μ' (mean) must be 1D for {dist_name}, got {mean_val.ndim}D.")
+ return None
+ if cov_val.ndim != 2:
+ logger.error(f"MvNormal 'Σ' (covariance) must be 2D for {dist_name}, got {cov_val.ndim}D.")
+ return None
+ if mean_val.shape[0] != cov_val.shape[0] or cov_val.shape[0] != cov_val.shape[1]:
+ logger.error(f"MvNormal mean/covariance dimension mismatch for {dist_name}: mean {mean_val.shape}, cov {cov_val.shape}")
+ return None
+ try:
+ return stats.multivariate_normal(mean=mean_val, cov=cov_val, allow_singular=True)
+ except ValueError as ve:
+ logger.error(f"Scipy MvNormal error for {dist_name} with mean {mean_val} and cov {cov_val}: {ve}")
+ return None
+ return None
+
+ def _convert_mv_tdist():
+ # Distributions.jl: MvTDist(df, μ, Σ) where Σ is scale matrix (not covariance)
+ # scipy.stats.multivariate_t(loc=μ, shape=Σ, df=df) where shape is scale matrix
+ df = _get_scalar_param("df") # degrees of freedom (from Distributions.jl MvTDist struct)
+ # Note: Distributions.jl MvTDist constructor is MvTDist(ν, μ, Σ), but fields might be df, loc, invscale, etc.
+ # Assuming direct field access for μ and Σ based on common parameter names.
+ # If Distributions.jl stores them as `loc` and `scale` for MvTDist, adjust "μ" and "Σ".
+ # For now, trying "μ" and "Σ" as per the original map, and "df" for degrees of freedom.
+ # If "df" is not found, try "ν" as it's common for T-distributions.
+ if df is None:
+ df = _get_scalar_param("ν")
+
+ mu_loc = _get_numpy_array_param("μ") # location vector
+ sigma_shape = _get_numpy_array_param("Σ") # scale matrix
+
+ if df is not None and mu_loc is not None and sigma_shape is not None:
+ if df <= 0:
+ logger.error(f"MvTDist 'df' (degrees of freedom) must be positive for {dist_name}, got {df}.")
+ return None
+ if mu_loc.ndim != 1:
+ logger.error(f"MvTDist 'μ' (loc) must be 1D for {dist_name}, got {mu_loc.ndim}D.")
+ return None
+ if sigma_shape.ndim != 2:
+ logger.error(f"MvTDist 'Σ' (shape) must be 2D for {dist_name}, got {sigma_shape.ndim}D.")
+ return None
+ if mu_loc.shape[0] != sigma_shape.shape[0] or sigma_shape.shape[0] != sigma_shape.shape[1]:
+ logger.error(f"MvTDist loc/shape dimension mismatch for {dist_name}: loc {mu_loc.shape}, shape {sigma_shape.shape}")
+ return None
+ try:
+ # allow_singular might not be an option for multivariate_t, check scipy docs if issues arise
+ return stats.multivariate_t(loc=mu_loc, shape=sigma_shape, df=df)
+ except ValueError as ve:
+ logger.error(f"Scipy MvTDist error for {dist_name}: {ve}")
+ return None
+ return None
+
+ def _convert_dirichlet():
+ alpha_conc = _get_numpy_array_param("alpha") # concentration parameters
+ # Distributions.jl uses 'α', let's try that first.
+ if alpha_conc is None:
+ alpha_conc = _get_numpy_array_param("α")
+
+ if alpha_conc is not None:
+ if alpha_conc.ndim != 1:
+ logger.error(f"Dirichlet 'α' (concentration) must be 1D for {dist_name}, got {alpha_conc.ndim}D.")
+ return None
+ if not np.all(alpha_conc > 0):
+ logger.error(f"Dirichlet 'α' elements must be positive for {dist_name}, got {alpha_conc}.")
+ return None
+ return stats.dirichlet(alpha=alpha_conc)
+ return None
+
+ def _convert_bernoulli():
+ p = _get_scalar_param("p")
+ if p is not None:
+ if not (0 <= p <= 1):
+ logger.error(f"Bernoulli 'p' must be in [0,1] for {dist_name}, got {p}.")
+ return None
+ return stats.bernoulli(p=p)
+ return None
+
+ def _convert_categorical():
+ # Distributions.jl Categorical has 'p' (probability vector)
+ # scipy.stats.multinomial(n=1, p=prob_vector)
+ p_probs = _get_numpy_array_param("p")
+ if p_probs is not None:
+ if p_probs.ndim != 1:
+ logger.error(f"Categorical 'p' (probabilities) must be 1D for {dist_name}, got {p_probs.ndim}D.")
+ return None
+ if not np.all(p_probs >= 0):
+ logger.error(f"Categorical 'p' elements must be non-negative for {dist_name}, got {p_probs}.")
+ return None
+ if not np.isclose(np.sum(p_probs), 1.0):
+ logger.warning(f"Categorical 'p' elements for {dist_name} do not sum to 1 (sum={np.sum(p_probs)}). Normalizing for scipy.")
+ p_probs = p_probs / np.sum(p_probs)
+ return stats.multinomial(n=1, p=p_probs)
+ return None
+
+ def _convert_multinomial():
+ n_trials = _get_int_param("n")
+ p_probs = _get_numpy_array_param("p")
+ if n_trials is not None and p_probs is not None:
+ if n_trials < 0: # typically n_trials >= 1 for multinomial, but scipy allows n=0
+ logger.error(f"Multinomial 'n' (trials) must be non-negative for {dist_name}, got {n_trials}.")
+ return None
+ if p_probs.ndim != 1:
+ logger.error(f"Multinomial 'p' (probabilities) must be 1D for {dist_name}, got {p_probs.ndim}D.")
+ return None
+ if not np.all(p_probs >= 0):
+ logger.error(f"Multinomial 'p' elements must be non-negative for {dist_name}, got {p_probs}.")
+ return None
+ if not np.isclose(np.sum(p_probs), 1.0):
+ logger.warning(f"Multinomial 'p' elements for {dist_name} do not sum to 1 (sum={np.sum(p_probs)}). Normalizing for scipy.")
+ p_probs = p_probs / np.sum(p_probs)
+ return stats.multinomial(n=n_trials, p=p_probs)
+ return None
+
+ # --- Conversion Map ---
+ conversion_map = {
+ "Normal": _convert_normal,
+ "Exponential": _convert_exponential,
+ "Uniform": _convert_uniform,
+ "Gamma": _convert_gamma,
+ "Beta": _convert_beta,
+ "Binomial": _convert_binomial,
+ "Poisson": _convert_poisson,
+ "Geometric": _convert_geometric,
+ "NegativeBinomial": _convert_negative_binomial,
+ "ChiSquared": _convert_chi_squared, # Julia: ChiSquared, Scipy: chi2
+ "Chisq": _convert_chi_squared, # Alias for ChiSquared
+ "TDist": _convert_tdist,
+ "Laplace": _convert_laplace,
+ "Cauchy": _convert_cauchy,
+ "LogNormal": _convert_log_normal,
+ "Weibull": _convert_weibull,
+ "Gumbel": _convert_gumbel,
+ "MvNormal": _convert_mv_normal,
+ "MvTDist": _convert_mv_tdist,
+ "Dirichlet": _convert_dirichlet,
+ "Bernoulli": _convert_bernoulli,
+ "Categorical": _convert_categorical,
+ "Multinomial": _convert_multinomial,
+ # Add InverseGamma from distributions.py if needed
+ "InverseGamma": lambda: ( # Distributions.jl: InverseGamma(shape, scale) = InverseGamma(α, θ)
+ alpha_ig := _get_scalar_param("α"), # shape
+ theta_ig := _get_scalar_param("θ"), # scale
+ (stats.invgamma(a=alpha_ig, scale=theta_ig)
+ if alpha_ig is not None and theta_ig is not None and alpha_ig > 0 and theta_ig > 0
+ else (logger.error(f"InverseGamma 'α' (shape) and 'θ' (scale) must be positive for {dist_name}, got α={alpha_ig}, θ={theta_ig}."), None)[1] if alpha_ig is not None and theta_ig is not None else None)
+ ),
+ # Add FDist from distributions.py if needed
+ "FDist": lambda: (
+ dfn := _get_scalar_param("ν1"),
+ dfd := _get_scalar_param("ν2"),
+ (stats.f(dfn=dfn, dfd=dfd)
+ if dfn is not None and dfd is not None and dfn > 0 and dfd > 0
+ else (logger.error(f"FDist 'ν1' and 'ν2' must be positive for {dist_name}, got ν1={dfn}, ν2={dfd}."), None)[1] if dfn is not None and dfd is not None else None)
+ )
+ }
+
+ if dist_name in conversion_map:
+ return conversion_map[dist_name]()
+ else:
+ logger.warning(f"Unsupported Julia distribution type for scipy conversion: {dist_name}")
+ return None
+
+ except Exception as e:
+ # Use dist_name_for_log which might be "UnknownDistribution" or the actual name if obtained
+ logger.error(f"Error converting Julia distribution ({dist_name_for_log}) to scipy: {e}", exc_info=True)
+ return None
@dataclass
@@ -1077,7 +1338,8 @@ class EstimatedParameters:
)
@classmethod
- def from_julia(cls, jl_estimatedparameters) -> Self:
+ def from_julia(cls, jl_estimatedparameters: Any) -> "EstimatedParameters":
+ """Converts a Julia EstimatedParameters object to its Python equivalent."""
return EstimatedParameters(
index=list(jl_estimatedparameters.index),
initialvalue=list(jl_estimatedparameters.initialvalue),
@@ -1175,7 +1437,8 @@ class Work:
)
@classmethod
- def from_julia(cls, jl_work) -> Self:
+ def from_julia(cls, jl_work) -> "Work":
+ """Converts a Julia Work object to its Python equivalent."""
work = Work()
work.analytical_steadystate_variables = list(
jl_work.analytical_steadystate_variables
@@ -1226,7 +1489,8 @@ class DynareSymbol:
)
@classmethod
- def from_julia(cls, jl_dynaresymbol) -> Self:
+ def from_julia(cls, jl_dynaresymbol: Any) -> "DynareSymbol":
+ """Converts a Julia DynareSymbol object to its Python equivalent."""
return DynareSymbol(
longname=jl_dynaresymbol.longname,
texname=jl_dynaresymbol.texname,
@@ -1267,8 +1531,9 @@ class Context:
)
@classmethod
- def from_julia(cls, jl_context) -> Self:
- context_py = Context(
+ def from_julia(cls, jl_context: Any) -> "Context":
+ """Converts a Julia Context object to its Python equivalent."""
+ return Context(
symboltable=symboltable_from_julia(jl_context.symboltable),
models=[Model.from_julia(jlm) for jlm in jl_context.models],
modfileinfo=ModFileInfo.from_julia(jl_context.modfileinfo),
@@ -1277,4 +1542,3 @@ class Context:
workspaces=None,
# workspaces=from_julia(jl_context.workspaces),
)
- return context_py
diff --git a/src/dynare/errors.py b/src/dynare/errors.py
index 2dfb6a2..a10ea84 100644
--- a/src/dynare/errors.py
+++ b/src/dynare/errors.py
@@ -1,6 +1,3 @@
-from typing import Self
-
-
class DynareError(Exception):
"""Exception raised for errors occurring during the execution of the dynare Julia command."""
@@ -12,6 +9,6 @@ class DynareError(Exception):
return f"DynareError: {self.message}"
@classmethod
- def from_julia_error(cls, julia_error) -> Self:
+ def from_julia_error(cls, julia_error) -> "DynareError":
message = f"JuliaError:\n{str(julia_error)}"
return cls(message)
diff --git a/src/dynare/indices.py b/src/dynare/indices.py
new file mode 100644
index 0000000..a4862b5
--- /dev/null
+++ b/src/dynare/indices.py
@@ -0,0 +1,123 @@
+from dataclasses import dataclass, field
+from typing import List, Dict, Tuple, Any
+
+@dataclass
+class PyIndices:
+ current: List[int] = field(default_factory=list)
+ forward: List[int] = field(default_factory=list)
+ purely_forward: List[int] = field(default_factory=list)
+ backward: List[int] = field(default_factory=list)
+ both: List[int] = field(default_factory=list)
+ non_backward: List[int] = field(default_factory=list)
+ static: List[int] = field(default_factory=list)
+ dynamic: List[int] = field(default_factory=list)
+ dynamic_current: List[int] = field(default_factory=list)
+ current_in_dynamic: List[int] = field(default_factory=list)
+ forward_in_dynamic: List[int] = field(default_factory=list)
+ backward_in_dynamic: List[int] = field(default_factory=list)
+ current_in_dynamic_jacobian: List[int] = field(default_factory=list)
+ current_in_static_jacobian: List[int] = field(default_factory=list)
+ exogenous: List[int] = field(default_factory=list)
+ n_endogenous: int = 0
+ D_columns: Dict[str, List[int]] = field(default_factory=dict)
+ E_columns: Dict[str, List[int]] = field(default_factory=dict)
+ UD_columns: List[int] = field(default_factory=list)
+ UE_columns: List[int] = field(default_factory=list)
+
+ @classmethod
+ def from_julia(cls, julia_indices: Any) -> "PyIndices":
+ if julia_indices is None:
+ return cls()
+
+ d_cols_data = {}
+ jl_D_columns = getattr(julia_indices, 'D_columns', None)
+ if jl_D_columns is not None:
+ if hasattr(jl_D_columns, 'D') and hasattr(jl_D_columns, 'jacobian'):
+ d_cols_data = {
+ "D": list(jl_D_columns.D) if jl_D_columns.D is not None else [],
+ "jacobian": list(jl_D_columns.jacobian) if jl_D_columns.jacobian is not None else [],
+ }
+ elif isinstance(jl_D_columns, tuple) and len(jl_D_columns) == 2: # Handles pre-converted tuple
+ d_cols_data = {
+ "D": list(jl_D_columns[0]) if jl_D_columns[0] is not None else [],
+ "jacobian": list(jl_D_columns[1]) if jl_D_columns[1] is not None else [],
+ }
+
+ e_cols_data = {}
+ jl_E_columns = getattr(julia_indices, 'E_columns', None)
+ if jl_E_columns is not None:
+ if hasattr(jl_E_columns, 'E') and hasattr(jl_E_columns, 'jacobian'):
+ e_cols_data = {
+ "E": list(jl_E_columns.E) if jl_E_columns.E is not None else [],
+ "jacobian": list(jl_E_columns.jacobian) if jl_E_columns.jacobian is not None else [],
+ }
+ elif isinstance(jl_E_columns, tuple) and len(jl_E_columns) == 2: # Handles pre-converted tuple
+ e_cols_data = {
+ "E": list(jl_E_columns[0]) if jl_E_columns[0] is not None else [],
+ "jacobian": list(jl_E_columns[1]) if jl_E_columns[1] is not None else [],
+ }
+
+ # Helper to safely convert attributes that are expected to be lists
+ def _to_list(attr_name: str) -> List[Any]:
+ attr_val = getattr(julia_indices, attr_name, [])
+ if attr_val is None:
+ return []
+ try:
+ return list(attr_val)
+ except TypeError: # If it's not iterable (e.g. an int like n_endogenous)
+ return []
+
+ return cls(
+ current=_to_list('current'),
+ forward=_to_list('forward'),
+ purely_forward=_to_list('purely_forward'),
+ backward=_to_list('backward'),
+ both=_to_list('both'),
+ non_backward=_to_list('non_backward'),
+ static=_to_list('static'),
+ dynamic=_to_list('dynamic'),
+ dynamic_current=_to_list('dynamic_current'),
+ current_in_dynamic=_to_list('current_in_dynamic'),
+ forward_in_dynamic=_to_list('forward_in_dynamic'),
+ backward_in_dynamic=_to_list('backward_in_dynamic'),
+ current_in_dynamic_jacobian=_to_list('current_in_dynamic_jacobian'),
+ current_in_static_jacobian=_to_list('current_in_static_jacobian'),
+ exogenous=_to_list('exogenous'),
+ n_endogenous=getattr(julia_indices, 'n_endogenous', 0),
+ D_columns=d_cols_data,
+ E_columns=e_cols_data,
+ UD_columns=_to_list('UD_columns'),
+ UE_columns=_to_list('UE_columns'),
+ )
+
+ def get_D_columns_as_tuple(self) -> Tuple[List[int], List[int]]:
+ return (self.D_columns.get("D", []), self.D_columns.get("jacobian", []))
+
+ def get_E_columns_as_tuple(self) -> Tuple[List[int], List[int]]:
+ return (self.E_columns.get("E", []), self.E_columns.get("jacobian", []))
+
+ def __repr__(self) -> str:
+ return (
+ f"PyIndices(\n"
+ f" current={self.current},\n"
+ f" forward={self.forward},\n"
+ f" purely_forward={self.purely_forward},\n"
+ f" backward={self.backward},\n"
+ f" both={self.both},\n"
+ f" non_backward={self.non_backward},\n"
+ f" static={self.static},\n"
+ f" dynamic={self.dynamic},\n"
+ f" dynamic_current={self.dynamic_current},\n"
+ f" current_in_dynamic={self.current_in_dynamic},\n"
+ f" forward_in_dynamic={self.forward_in_dynamic},\n"
+ f" backward_in_dynamic={self.backward_in_dynamic},\n"
+ f" current_in_dynamic_jacobian={self.current_in_dynamic_jacobian},\n"
+ f" current_in_static_jacobian={self.current_in_static_jacobian},\n"
+ f" exogenous={self.exogenous},\n"
+ f" n_endogenous={self.n_endogenous},\n"
+ f" D_columns={self.D_columns},\n"
+ f" E_columns={self.E_columns},\n"
+ f" UD_columns={self.UD_columns},\n"
+ f" UE_columns={self.UE_columns}\n"
+ f")"
+ )
diff --git a/src/dynare/juliapkg.json b/src/dynare/juliapkg.json
index 6464d88..9556035 100644
--- a/src/dynare/juliapkg.json
+++ b/src/dynare/juliapkg.json
@@ -1,6 +1,10 @@
{
"julia": "1.11",
"packages": {
+ "TransformVariables": {
+ "uuid": "84d833dd-6860-57f9-a1a7-6da5db126cff",
+ "version": "=0.8.10"
+ },
"DataFrames": {
"uuid": "a93c6f00-e57d-5684-b7b6-d8193f3e46c0",
"version": "1.7.0"
diff --git a/src/dynare/workspace.py b/src/dynare/workspace.py
index 3ed2942..2ab1c49 100644
--- a/src/dynare/workspace.py
+++ b/src/dynare/workspace.py
@@ -1,61 +1,7 @@
-from dataclasses import dataclass, field
-from typing import List, Any, Dict, Tuple, Union, Self, NamedTuple
+from dataclasses import dataclass
+from typing import List, Union, NamedTuple, Any
import numpy as np
-
-
-@dataclass
-class PyIndices:
- current: List[int]
- forward: List[int]
- purely_forward: List[int]
- backward: List[int]
- both: List[int]
- non_backward: List[int]
- static: List[int]
- dynamic: List[int]
- dynamic_current: List[int]
- current_in_dynamic: List[int]
- forward_in_dynamic: List[int]
- backward_in_dynamic: List[int]
- current_in_dynamic_jacobian: List[int]
- current_in_static_jacobian: List[int]
- exogenous: List[int]
- n_endogenous: int
- D_columns: NamedTuple
- E_columns: NamedTuple
- UD_columns: List[int]
- UE_columns: List[int]
-
- @classmethod
- def from_julia(cls, julia_indices) -> Self:
- return cls(
- current=list(julia_indices.current),
- forward=list(julia_indices.forward),
- purely_forward=list(julia_indices.purely_forward),
- backward=list(julia_indices.backward),
- both=list(julia_indices.both),
- non_backward=list(julia_indices.non_backward),
- static=list(julia_indices.static),
- dynamic=list(julia_indices.dynamic),
- dynamic_current=list(julia_indices.dynamic_current),
- current_in_dynamic=list(julia_indices.current_in_dynamic),
- forward_in_dynamic=list(julia_indices.forward_in_dynamic),
- backward_in_dynamic=list(julia_indices.backward_in_dynamic),
- current_in_dynamic_jacobian=list(julia_indices.current_in_dynamic_jacobian),
- current_in_static_jacobian=list(julia_indices.current_in_static_jacobian),
- exogenous=list(julia_indices.exogenous),
- n_endogenous=julia_indices.n_endogenous,
- D_columns={
- "D": list(julia_indices.D_columns.D),
- "jacobian": list(julia_indices.D_columns.jacobian),
- },
- E_columns={
- "E": list(julia_indices.E_columns.E),
- "jacobian": list(julia_indices.E_columns.jacobian),
- },
- UD_columns=list(julia_indices.UD_columns),
- UE_columns=list(julia_indices.UE_columns),
- )
+from .indices import PyIndices # Added import
@dataclass
@@ -64,7 +10,8 @@ class PyCyclicReductionOptions:
tol: float = 1e-8
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyCyclicReductionOptions":
+ """Converts a Julia CyclicReductionOptions object to its Python equivalent."""
return cls(maxiter=julia_obj.maxiter, tol=julia_obj.tol)
@@ -73,7 +20,8 @@ class PyGeneralizedSchurOptions:
criterium: float = 1.0 + 1e-6
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyGeneralizedSchurOptions":
+ """Converts a Julia GeneralizedSchurOptions object to its Python equivalent."""
return cls(criterium=julia_obj.criterium)
@@ -83,7 +31,8 @@ class PyLinearRationalExpectationsOptions:
generalized_schur: PyGeneralizedSchurOptions = PyGeneralizedSchurOptions()
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyLinearRationalExpectationsOptions":
+ """Converts a Julia LinearRationalExpectationsOptions object to its Python equivalent."""
return cls(
cyclic_reduction=PyCyclicReductionOptions.from_julia(
julia_obj.cyclic_reduction
@@ -105,7 +54,8 @@ class PySchurWs:
eigen_values: List[complex]
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PySchurWs":
+ """Converts a Julia SchurWs object to its Python equivalent PySchurWs."""
return cls(
work=list(julia_obj.work),
wr=list(julia_obj.wr),
@@ -123,7 +73,8 @@ class PyLSEWs:
X: List[complex]
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyLSEWs":
+ """Converts a Julia LSEWs object to its Python equivalent PyLSEWs."""
return cls(work=julia_obj.work.to_numpy(), X=julia_obj.X.to_numpy())
@@ -140,7 +91,8 @@ class PyEigenWs:
rcondv: List[float]
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyEigenWs":
+ """Converts a Julia EigenWs object to its Python equivalent PyEigenWs."""
return cls(
work=list(julia_obj.work),
rwork=list(julia_obj.rwork),
@@ -164,7 +116,8 @@ class PyHermitianEigenWs:
isuppz: List
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyHermitianEigenWs":
+ """Converts a Julia HermitianEigenWs object to its Python equivalent PyHermitianEigenWs."""
return cls(
work=list(julia_obj.work),
rwork=list(julia_obj.rwork),
@@ -177,10 +130,11 @@ class PyHermitianEigenWs:
@dataclass
class PyLUWs:
- ipiv: List(int)
+ ipiv: List[int]
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyLUWs":
+ """Converts a Julia LUWs object to its Python equivalent PyLUWs."""
return cls(ipiv=julia_obj.ipiv.to_numpy())
@@ -190,7 +144,8 @@ class PyQRWs:
τ: np.ndarray
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyQRWs":
+ """Converts a Julia QRWs object to its Python equivalent PyQRWs."""
return cls(work=julia_obj.work.to_numpy(), τ=julia_obj.τ.to_numpy())
@@ -200,7 +155,8 @@ class PyQRWYWs:
T: np.ndarray
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyQRWYWs":
+ """Converts a Julia QRWYWs object to its Python equivalent PyQRWYWs."""
return cls(work=julia_obj.work.to_numpy(), T=julia_obj.T.to_numpy())
@@ -212,7 +168,8 @@ class PyQRPivotedWs:
jpvt: np.ndarray
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyQRPivotedWs":
+ """Converts a Julia QRPivotedWs object to its Python equivalent PyQRPivotedWs."""
return cls(
work=julia_obj.work.to_numpy(),
rwork=julia_obj.rwork.to_numpy(),
@@ -227,7 +184,8 @@ class PyQROrmWs:
τ: np.ndarray
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyQROrmWs":
+ """Converts a Julia QROrmWs object to its Python equivalent PyQROrmWs."""
return cls(work=julia_obj.work.to_numpy(), τ=julia_obj.τ.to_numpy())
@@ -247,7 +205,8 @@ class PyLyapdWs:
linsolve_ws2: "PyLUWs"
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyLyapdWs":
+ """Converts a Julia LyapdWs object to its Python equivalent PyLyapdWs."""
return cls(
AA=julia_obj.AA.to_numpy(),
AAtemp=julia_obj.AAtemp.to_numpy(),
@@ -279,7 +238,8 @@ class PyNonstationaryVarianceWs:
nonstate_stationary_variables: List[bool]
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyNonstationaryVarianceWs":
+ """Converts a Julia NonstationaryVarianceWs object to its Python equivalent."""
return cls(
rΣ_s_s=julia_obj.rΣ_s_s.to_numpy(),
rA1=julia_obj.rA1.to_numpy(),
@@ -310,7 +270,8 @@ class PyVarianceWs:
lyapd_ws: PyLyapdWs
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyVarianceWs":
+ """Converts a Julia VarianceWs object to its Python equivalent PyVarianceWs."""
return cls(
B1S=julia_obj.B1S.to_numpy(),
B1SB1=julia_obj.B1SB1.to_numpy(),
@@ -335,7 +296,8 @@ class PyBunchKaufmanWs:
ipiv: np.ndarray
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyBunchKaufmanWs":
+ """Converts a Julia BunchKaufmanWs object to its Python equivalent."""
return cls(work=julia_obj.work.to_numpy(), ipiv=julia_obj.ipiv.to_numpy())
@@ -345,7 +307,8 @@ class PyCholeskyPivotedWs:
piv: np.ndarray
@classmethod
- def from_julia(cls, julia_obj):
+ def from_julia(cls, julia_obj: Any) -> "PyCholeskyPivotedWs":
+ """Converts a Julia CholeskyPivotedWs object to its Python equivalent."""
return cls(work=julia_obj.work.to_numpy(), piv=julia_obj.piv.to_numpy())
@@ -360,7 +323,17 @@ class PyGsSolverWs:
schurws: PyGeneralizedSchurWs
@classmethod
- def from_julia(cls, d, n1):
+ def from_julia(cls, d: Any, n1: int) -> "PyGsSolverWs":
+ """
+ Converts Julia GsSolverWs components to Python equivalent PyGsSolverWs.
+
+ Args:
+ d: Julia object, typically a matrix, used for schurws.
+ n1: Integer dimension.
+
+ Returns:
+ An instance of PyGsSolverWs.
+ """
n = d.shape[0]
n2 = n - n1
tmp1 = np.empty((n1, n1)) # Create similar matrix for tmp1
@@ -384,7 +357,8 @@ class PyLinearGsSolverWs:
e: List[float]
@classmethod
- def from_julia(cls, julia_solver_ws) -> "PyLinearGsSolverWs":
+ def from_julia(cls, julia_solver_ws: Any) -> "PyLinearGsSolverWs":
+ """Converts a Julia LinearGsSolverWs object to its Python equivalent."""
return cls(
solver_ws=PyGsSolverWs.from_julia(
julia_solver_ws.solver_ws
@@ -405,7 +379,8 @@ class PyLinearCyclicReductionWs:
x: List[float]
@classmethod
- def from_julia(cls, julia_cyclic_ws) -> "PyLinearCyclicReductionWs":
+ def from_julia(cls, julia_cyclic_ws: Any) -> "PyLinearCyclicReductionWs":
+ """Converts a Julia LinearCyclicReductionWs object to its Python equivalent."""
return cls(
solver_ws=PyCyclicReductionWs.from_julia(
julia_cyclic_ws.solver_ws
@@ -440,7 +415,8 @@ class PyLinearRationalExpectationsWs:
AGplusB_linsolve_ws: "PyLUWs"
@classmethod
- def from_julia(cls, julia_lre_ws) -> Self:
+ def from_julia(cls, julia_lre_ws: Any) -> "PyLinearRationalExpectationsWs":
+ """Converts a Julia LinearRationalExpectationsWs object to its Python equivalent."""
solver_ws = (
PyLinearGsSolverWs.from_julia(julia_lre_ws.solver_ws)
if type(julia_lre_ws.solver_ws) == "LinearGsSolverWs"
--
GitLab
From d9ef60fab005a5fc89cadc5a36dd688df7208f79 Mon Sep 17 00:00:00 2001
From: Daniel Sali <daniel.sali@alphacruncher.com>
Date: Tue, 3 Jun 2025 09:26:35 +0200
Subject: [PATCH 2/6] Add InverseGamma1 distribution support, black formatting
---
src/dynare/__init__.py | 7 +-
src/dynare/dynare.py | 18 +-
src/dynare/dynare_context.py | 407 +++++++++++++++++++++++++----------
src/dynare/indices.py | 81 ++++---
4 files changed, 358 insertions(+), 155 deletions(-)
diff --git a/src/dynare/__init__.py b/src/dynare/__init__.py
index 81b1582..d57010b 100644
--- a/src/dynare/__init__.py
+++ b/src/dynare/__init__.py
@@ -7,7 +7,10 @@ from .dynare import dynare
logger = logging.getLogger("DynarePython")
-def setup_logging(level_str: str = "INFO", format_string: str = "%(asctime)s - %(name)s - %(levelname)s - %(message)s"):
+def setup_logging(
+ level_str: str = "INFO",
+ format_string: str = "%(asctime)s - %(name)s - %(levelname)s - %(message)s",
+):
"""
Set up logging for the DynarePython library.
@@ -62,4 +65,4 @@ def _run_once():
_run_once()
-__all__ = ["dynare", "setup_logging"]
\ No newline at end of file
+__all__ = ["dynare", "setup_logging"]
diff --git a/src/dynare/dynare.py b/src/dynare/dynare.py
index aa73fb9..a097bab 100644
--- a/src/dynare/dynare.py
+++ b/src/dynare/dynare.py
@@ -14,15 +14,19 @@ def dynare(model: str | Path, *dynare_options: str) -> Context:
model = Path(model)
jl.seval("using Serialization")
jl.seval("using Dynare")
-
+
resolved_path_str = str(model.resolve())
- escaped_path_for_julia_literal = resolved_path_str.replace('\\', '\\\\').replace('"', '\\"')
+ escaped_path_for_julia_literal = resolved_path_str.replace(
+ "\\", "\\\\"
+ ).replace('"', '\\"')
julia_options_str_parts = []
for opt_str in dynare_options:
- escaped_opt_for_julia_literal = opt_str.replace('\\', '\\\\').replace('"', '\\"')
+ escaped_opt_for_julia_literal = opt_str.replace("\\", "\\\\").replace(
+ '"', '\\"'
+ )
julia_options_str_parts.append(f'"{escaped_opt_for_julia_literal}"')
-
+
julia_options_cmd_part = " ".join(julia_options_str_parts)
jl.seval(
@@ -185,14 +189,14 @@ def dynare(model: str | Path, *dynare_options: str) -> Context:
end
"""
)
-
- julia_command = f'''ctx = @dynare "{escaped_path_for_julia_literal}" {julia_options_cmd_part};
+
+ julia_command = f"""ctx = @dynare "{escaped_path_for_julia_literal}" {julia_options_cmd_part};
if !(ctx isa Dynare.Context)
throw(error("Failed to produce a Dynare context."))
else
convert_to_pycontext(ctx)
end
- '''
+ """
context = jl.seval(julia_command)
return Context.from_julia(context)
except JuliaError as e:
diff --git a/src/dynare/dynare_context.py b/src/dynare/dynare_context.py
index ea24492..aa2bd00 100644
--- a/src/dynare/dynare_context.py
+++ b/src/dynare/dynare_context.py
@@ -7,6 +7,7 @@ from typing import List, Any, Dict, Tuple, Union
import numpy as np
import pandas as pd
import scipy.stats as stats
+from scipy.special import gamma as _gamma_func # for InverseGamma1
from juliacall import Main as jl # Added import
from .indices import PyIndices
@@ -338,9 +339,11 @@ class Model:
),
parameter_nbr=get_member(jl_model, "parameter_nbr"),
original_endogenous_nbr=get_member(jl_model, "original_endogenous_nbr"),
- lead_lag_incidence=get_member(jl_model, "lead_lag_incidence").to_numpy()
- if get_member(jl_model, "lead_lag_incidence")
- else None,
+ lead_lag_incidence=(
+ get_member(jl_model, "lead_lag_incidence").to_numpy()
+ if get_member(jl_model, "lead_lag_incidence")
+ else None
+ ),
n_static=get_member(jl_model, "n_static"),
n_fwrd=get_member(jl_model, "n_fwrd"),
n_bkwrd=get_member(jl_model, "n_bkwrd"),
@@ -350,9 +353,11 @@ class Model:
DErows2=list(get_member(jl_model, "DErows2", [])),
n_dyn=get_member(jl_model, "n_dyn"),
i_static=list(get_member(jl_model, "i_static", [])),
- i_dyn=get_member(jl_model, "i_dyn").to_numpy()
- if get_member(jl_model, "i_dyn")
- else None,
+ i_dyn=(
+ get_member(jl_model, "i_dyn").to_numpy()
+ if get_member(jl_model, "i_dyn")
+ else None
+ ),
i_bkwrd=list(get_member(jl_model, "i_bkwrd", [])),
i_bkwrd_b=list(get_member(jl_model, "i_bkwrd_b", [])),
i_bkwrd_ns=list(get_member(jl_model, "i_bkwrd_ns", [])),
@@ -396,9 +401,11 @@ class Model:
serially_correlated_exogenous=list(
get_member(jl_model, "serially_correlated_exogenous", [])
),
- Sigma_e=get_member(jl_model, "Sigma_e").to_numpy()
- if get_member(jl_model, "Sigma_e")
- else None,
+ Sigma_e=(
+ get_member(jl_model, "Sigma_e").to_numpy()
+ if get_member(jl_model, "Sigma_e")
+ else None
+ ),
maximum_endo_lag=get_member(jl_model, "maximum_endo_lag"),
maximum_endo_lead=get_member(jl_model, "maximum_endo_lead"),
maximum_exo_lag=get_member(jl_model, "maximum_exo_lag"),
@@ -623,7 +630,9 @@ class LinearRationalExpectationsResults:
)
@classmethod
- def from_julia(cls, jl_linearrationalexpectationsresults: Any) -> "LinearRationalExpectationsResults":
+ def from_julia(
+ cls, jl_linearrationalexpectationsresults: Any
+ ) -> "LinearRationalExpectationsResults":
"""Converts a Julia LinearRationalExpectationsResults object to its Python equivalent."""
return LinearRationalExpectationsResults(
eigenvalues=jl_linearrationalexpectationsresults.eigenvalues.to_numpy(),
@@ -683,7 +692,9 @@ class PyTasmanianSG:
def from_julia(cls, jl_tasmaniansg: Any) -> "PyTasmanianSG":
"""Converts a Julia TasmanianSG object to its Python equivalent PyTasmanianSG."""
if not TASlib:
- raise ValueError("It looks like Tasmanian is not installed. Please install from https://github.com/ORNL/TASMANIAN and set LIBTASMANIAN_PATH")
+ raise ValueError(
+ "It looks like Tasmanian is not installed. Please install from https://github.com/ORNL/TASMANIAN and set LIBTASMANIAN_PATH"
+ )
c_func = getattr(TASlib, "tsgConstructTasmanianSparseGrid")
c_func.restype = ctypes.POINTER(ctypes.c_void_p)
pGrid = c_func()
@@ -699,23 +710,27 @@ class PyTasmanianSG:
)
def __del__(self):
- if hasattr(self, 'pGrid') and self.pGrid and TASlib:
+ if hasattr(self, "pGrid") and self.pGrid and TASlib:
try:
# Attempt to get the destructor function from the loaded TASlib
destructor_func = getattr(TASlib, "tsgDestructTasmanianSparseGrid")
destructor_func.argtypes = [ctypes.POINTER(ctypes.c_void_p)]
destructor_func.restype = None # Typically void for destructors
-
+
destructor_func(self.pGrid)
- logger.debug(f"Tasmanian sparse grid {self.pGrid} destructed successfully.")
+ logger.debug(
+ f"Tasmanian sparse grid {self.pGrid} destructed successfully."
+ )
except AttributeError:
# This might happen if TASlib is loaded but doesn't have the expected destructor
logger.error(
"Tasmanian library is loaded but 'tsgDestructTasmanianSparseGrid' function not found. Cannot free pGrid."
)
except Exception as e:
- logger.error(f"Error during destruction of Tasmanian sparse grid {self.pGrid}: {e}")
- elif hasattr(self, 'pGrid') and self.pGrid and not TASlib:
+ logger.error(
+ f"Error during destruction of Tasmanian sparse grid {self.pGrid}: {e}"
+ )
+ elif hasattr(self, "pGrid") and self.pGrid and not TASlib:
logger.warning(
f"Tasmanian sparse grid {self.pGrid} might not be freed because TASlib is not loaded (LIBTASMANIAN_PATH not set or invalid)."
)
@@ -843,7 +858,9 @@ class ModelResults:
forecast: List[
pd.DataFrame
] # Assuming AxisArrayTable can be represented as a DataFrame
- initial_smoother: pd.DataFrame # Assuming AxisArrayTable can be represented as a DataFrame
+ initial_smoother: (
+ pd.DataFrame
+ ) # Assuming AxisArrayTable can be represented as a DataFrame
linearrationalexpectations: LinearRationalExpectationsResults
simulations: List[Simulation]
smoother: pd.DataFrame # Assuming AxisArrayTable can be represented as a DataFrame
@@ -925,8 +942,36 @@ def distribution_from_julia(
]:
"""
Converts a Julia Distributions.jl distribution object to its equivalent scipy.stats object.
- Enhanced for robustness, safer Julia interaction, and detailed validation.
"""
+
+ class _InvGamma1(stats.rv_continuous):
+ """
+ The *inverse Gamma distribution* with shape parameter `α` and scale `θ` has probability
+ density function
+
+ ```math
+ f(x; \\alpha, \\theta) = 2\\frac{\\theta^\\alpha x^{-(2\\alpha + 1)}}{\\Gamma(\\alpha)}
+ e^{-\\frac{\\theta}{x^2}}, \\quad x > 0
+ ```
+ """
+
+ def _pdf(self, x, a):
+ return (
+ 2
+ * self._scale**a
+ / _gamma_func(a)
+ * x ** (-2 * a - 1)
+ * np.exp(-self._scale / (x**2))
+ )
+
+ def _rvs(self, a, size=None, random_state=None):
+ y = stats.invgamma(a=a, scale=self._scale).rvs(
+ size=size, random_state=random_state
+ )
+ return np.sqrt(y)
+
+ _inv_gamma1 = _InvGamma1(name="inv_gamma1")
+
if jl_distribution is None:
return None
@@ -934,7 +979,7 @@ def distribution_from_julia(
try:
julia_type_obj = jl.typeof(jl_distribution)
dist_name = jl.nameof(julia_type_obj).__str__()
- dist_name_for_log = dist_name # Update for logging once name is known
+ dist_name_for_log = dist_name
# --- Helper functions for parameter extraction ---
def _get_scalar_param(param_name_str: str) -> Union[float, None]:
@@ -942,7 +987,9 @@ def distribution_from_julia(
jl_val = jl.getproperty(jl_distribution, jl.Symbol(param_name_str))
return jl.pyconvert(float, jl_val)
except Exception:
- logger.debug(f"Scalar parameter '{param_name_str}' not found or not convertible for {dist_name}.")
+ logger.debug(
+ f"Scalar parameter '{param_name_str}' not found or not convertible for {dist_name}."
+ )
return None
def _get_int_param(param_name_str: str) -> Union[int, None]:
@@ -950,7 +997,9 @@ def distribution_from_julia(
jl_val = jl.getproperty(jl_distribution, jl.Symbol(param_name_str))
return jl.pyconvert(int, jl_val)
except Exception:
- logger.debug(f"Integer parameter '{param_name_str}' not found or not convertible for {dist_name}.")
+ logger.debug(
+ f"Integer parameter '{param_name_str}' not found or not convertible for {dist_name}."
+ )
return None
def _get_numpy_array_param(param_name_str: str) -> Union[np.ndarray, None]:
@@ -958,7 +1007,9 @@ def distribution_from_julia(
jl_val = jl.getproperty(jl_distribution, jl.Symbol(param_name_str))
return jl.pyconvert(np.ndarray, jl_val)
except Exception:
- logger.debug(f"NumPy array parameter '{param_name_str}' not found or not convertible for {dist_name}.")
+ logger.debug(
+ f"NumPy array parameter '{param_name_str}' not found or not convertible for {dist_name}."
+ )
return None
# --- Individual conversion functions for each distribution type ---
@@ -967,18 +1018,22 @@ def distribution_from_julia(
sigma = _get_scalar_param("σ")
if mu is not None and sigma is not None:
if sigma <= 0:
- logger.error(f"Normal 'σ' must be positive for {dist_name}, got {sigma}.")
+ logger.error(
+ f"Normal 'σ' must be positive for {dist_name}, got {sigma}."
+ )
return None
return stats.norm(loc=mu, scale=sigma)
return None
def _convert_exponential():
- theta = _get_scalar_param("θ") # rate in Distributions.jl
+ theta = _get_scalar_param("θ") # rate in Distributions.jl
if theta is not None:
if theta <= 0:
- logger.error(f"Exponential 'θ' (rate) must be positive for {dist_name}, got {theta}.")
+ logger.error(
+ f"Exponential 'θ' (rate) must be positive for {dist_name}, got {theta}."
+ )
return None
- return stats.expon(scale=1/theta) # scipy.expon uses scale = 1/rate
+ return stats.expon(scale=1 / theta) # scipy.expon uses scale = 1/rate
return None
def _convert_uniform():
@@ -986,27 +1041,35 @@ def distribution_from_julia(
b = _get_scalar_param("b")
if a is not None and b is not None:
if b < a:
- logger.error(f"Uniform 'b' must be >= 'a' for {dist_name}, got a={a}, b={b}.")
+ logger.error(
+ f"Uniform 'b' must be >= 'a' for {dist_name}, got a={a}, b={b}."
+ )
return None
return stats.uniform(loc=a, scale=b - a)
return None
- def _convert_gamma(): # Distributions.jl: Gamma(shape, scale) = Gamma(α, θ)
- alpha = _get_scalar_param("α") # shape
- theta = _get_scalar_param("θ") # scale
+ def _convert_gamma(): # Distributions.jl: Gamma(shape, scale) = Gamma(α, θ)
+ alpha = _get_scalar_param("α") # shape
+ theta = _get_scalar_param("θ") # scale
if alpha is not None and theta is not None:
if alpha <= 0 or theta <= 0:
- logger.error(f"Gamma 'α' (shape) and 'θ' (scale) must be positive for {dist_name}, got α={alpha}, θ={theta}.")
+ logger.error(
+ f"Gamma 'α' (shape) and 'θ' (scale) must be positive for {dist_name}, got α={alpha}, θ={theta}."
+ )
return None
- return stats.gamma(a=alpha, scale=theta) # scipy.stats.gamma(a=shape, scale=scale)
+ return stats.gamma(
+ a=alpha, scale=theta
+ ) # scipy.stats.gamma(a=shape, scale=scale)
return None
-
+
def _convert_beta():
alpha = _get_scalar_param("α")
beta_p = _get_scalar_param("β")
if alpha is not None and beta_p is not None:
if alpha <= 0 or beta_p <= 0:
- logger.error(f"Beta 'α' and 'β' must be positive for {dist_name}, got α={alpha}, β={beta_p}.")
+ logger.error(
+ f"Beta 'α' and 'β' must be positive for {dist_name}, got α={alpha}, β={beta_p}."
+ )
return None
return stats.beta(a=alpha, b=beta_p)
return None
@@ -1016,10 +1079,14 @@ def distribution_from_julia(
p = _get_scalar_param("p")
if n is not None and p is not None:
if n < 0:
- logger.error(f"Binomial 'n' must be non-negative for {dist_name}, got n={n}.")
+ logger.error(
+ f"Binomial 'n' must be non-negative for {dist_name}, got n={n}."
+ )
return None
if not (0 <= p <= 1):
- logger.error(f"Binomial 'p' must be in [0,1] for {dist_name}, got p={p}.")
+ logger.error(
+ f"Binomial 'p' must be in [0,1] for {dist_name}, got p={p}."
+ )
return None
return stats.binom(n=n, p=p)
return None
@@ -1027,8 +1094,10 @@ def distribution_from_julia(
def _convert_poisson():
lambda_ = _get_scalar_param("λ")
if lambda_ is not None:
- if lambda_ < 0: # scipy allows lambda_ == 0
- logger.error(f"Poisson 'λ' must be non-negative for {dist_name}, got {lambda_}.")
+ if lambda_ < 0: # scipy allows lambda_ == 0
+ logger.error(
+ f"Poisson 'λ' must be non-negative for {dist_name}, got {lambda_}."
+ )
return None
return stats.poisson(mu=lambda_)
return None
@@ -1037,116 +1106,153 @@ def distribution_from_julia(
p = _get_scalar_param("p")
if p is not None:
if not (0 < p <= 1):
- logger.error(f"Geometric 'p' must be in (0,1] for {dist_name}, got {p}.")
+ logger.error(
+ f"Geometric 'p' must be in (0,1] for {dist_name}, got {p}."
+ )
return None
return stats.geom(p=p)
return None
def _convert_negative_binomial():
- r = _get_int_param("r") # number of successes
- p = _get_scalar_param("p") # probability of success
+ r = _get_int_param("r") # number of successes
+ p = _get_scalar_param("p") # probability of success
if r is not None and p is not None:
if r <= 0:
- logger.error(f"NegativeBinomial 'r' (successes) must be positive for {dist_name}, got {r}.")
+ logger.error(
+ f"NegativeBinomial 'r' (successes) must be positive for {dist_name}, got {r}."
+ )
return None
- if not (0 < p <= 1): # p can be 1
- logger.error(f"NegativeBinomial 'p' must be in (0,1] for {dist_name}, got {p}.")
+ if not (0 < p <= 1): # p can be 1
+ logger.error(
+ f"NegativeBinomial 'p' must be in (0,1] for {dist_name}, got {p}."
+ )
return None
return stats.nbinom(n=r, p=p)
return None
- def _convert_chi_squared(): # Julia: ChiSquared, Scipy: chi2
- df = _get_scalar_param("ν") # degrees of freedom
+ def _convert_chi_squared(): # Julia: ChiSquared, Scipy: chi2
+ df = _get_scalar_param("ν") # degrees of freedom
if df is not None:
if df <= 0:
- logger.error(f"ChiSquared 'ν' (df) must be positive for {dist_name}, got {df}.")
+ logger.error(
+ f"ChiSquared 'ν' (df) must be positive for {dist_name}, got {df}."
+ )
return None
return stats.chi2(df=df)
return None
-
+
def _convert_tdist():
- df = _get_scalar_param("ν") # degrees of freedom
+ df = _get_scalar_param("ν") # degrees of freedom
if df is not None:
if df <= 0:
- logger.error(f"TDist 'ν' (df) must be positive for {dist_name}, got {df}.")
+ logger.error(
+ f"TDist 'ν' (df) must be positive for {dist_name}, got {df}."
+ )
return None
return stats.t(df=df)
return None
def _convert_laplace():
- mu = _get_scalar_param("μ") # location
- beta_scale = _get_scalar_param("β") # scale
+ mu = _get_scalar_param("μ") # location
+ beta_scale = _get_scalar_param("β") # scale
if mu is not None and beta_scale is not None:
if beta_scale <= 0:
- logger.error(f"Laplace 'β' (scale) must be positive for {dist_name}, got {beta_scale}.")
+ logger.error(
+ f"Laplace 'β' (scale) must be positive for {dist_name}, got {beta_scale}."
+ )
return None
return stats.laplace(loc=mu, scale=beta_scale)
return None
def _convert_cauchy():
- x0 = _get_scalar_param("x0") # location
- gamma_scale = _get_scalar_param("γ") # scale
+ x0 = _get_scalar_param("x0") # location
+ gamma_scale = _get_scalar_param("γ") # scale
if x0 is not None and gamma_scale is not None:
if gamma_scale <= 0:
- logger.error(f"Cauchy 'γ' (scale) must be positive for {dist_name}, got {gamma_scale}.")
+ logger.error(
+ f"Cauchy 'γ' (scale) must be positive for {dist_name}, got {gamma_scale}."
+ )
return None
return stats.cauchy(loc=x0, scale=gamma_scale)
return None
def _convert_log_normal():
- mu_log = _get_scalar_param("μ") # log-mean
- sigma_log = _get_scalar_param("σ") # log-std
+ mu_log = _get_scalar_param("μ") # log-mean
+ sigma_log = _get_scalar_param("σ") # log-std
if mu_log is not None and sigma_log is not None:
if sigma_log <= 0:
- logger.error(f"LogNormal 'σ' (log-std) must be positive for {dist_name}, got {sigma_log}.")
+ logger.error(
+ f"LogNormal 'σ' (log-std) must be positive for {dist_name}, got {sigma_log}."
+ )
return None
- return stats.lognorm(s=sigma_log, scale=np.exp(mu_log)) # s is shape (log-std), scale is exp(log-mean)
+ return stats.lognorm(
+ s=sigma_log, scale=np.exp(mu_log)
+ ) # s is shape (log-std), scale is exp(log-mean)
return None
def _convert_weibull():
- k_shape = _get_scalar_param("k") # shape
- lambda_scale = _get_scalar_param("λ") # scale
+ k_shape = _get_scalar_param("k") # shape
+ lambda_scale = _get_scalar_param("λ") # scale
if k_shape is not None and lambda_scale is not None:
if k_shape <= 0 or lambda_scale <= 0:
- logger.error(f"Weibull 'k' (shape) and 'λ' (scale) must be positive for {dist_name}, got k={k_shape}, λ={lambda_scale}.")
+ logger.error(
+ f"Weibull 'k' (shape) and 'λ' (scale) must be positive for {dist_name}, got k={k_shape}, λ={lambda_scale}."
+ )
return None
return stats.weibull_min(c=k_shape, scale=lambda_scale)
return None
def _convert_gumbel():
- mu_loc = _get_scalar_param("μ") # location
- beta_scale = _get_scalar_param("β") # scale
+ mu_loc = _get_scalar_param("μ") # location
+ beta_scale = _get_scalar_param("β") # scale
if mu_loc is not None and beta_scale is not None:
if beta_scale <= 0:
- logger.error(f"Gumbel 'β' (scale) must be positive for {dist_name}, got {beta_scale}.")
+ logger.error(
+ f"Gumbel 'β' (scale) must be positive for {dist_name}, got {beta_scale}."
+ )
return None
return stats.gumbel_r(loc=mu_loc, scale=beta_scale)
return None
-
+
def _convert_mv_normal():
mean_val = _get_numpy_array_param("μ")
cov_val = _get_numpy_array_param("Σ")
if mean_val is not None and cov_val is not None:
if mean_val.ndim != 1:
- logger.error(f"MvNormal 'μ' (mean) must be 1D for {dist_name}, got {mean_val.ndim}D.")
+ logger.error(
+ f"MvNormal 'μ' (mean) must be 1D for {dist_name}, got {mean_val.ndim}D."
+ )
return None
if cov_val.ndim != 2:
- logger.error(f"MvNormal 'Σ' (covariance) must be 2D for {dist_name}, got {cov_val.ndim}D.")
+ logger.error(
+ f"MvNormal 'Σ' (covariance) must be 2D for {dist_name}, got {cov_val.ndim}D."
+ )
return None
- if mean_val.shape[0] != cov_val.shape[0] or cov_val.shape[0] != cov_val.shape[1]:
- logger.error(f"MvNormal mean/covariance dimension mismatch for {dist_name}: mean {mean_val.shape}, cov {cov_val.shape}")
+ if (
+ mean_val.shape[0] != cov_val.shape[0]
+ or cov_val.shape[0] != cov_val.shape[1]
+ ):
+ logger.error(
+ f"MvNormal mean/covariance dimension mismatch for {dist_name}: mean {mean_val.shape}, cov {cov_val.shape}"
+ )
return None
try:
- return stats.multivariate_normal(mean=mean_val, cov=cov_val, allow_singular=True)
+ return stats.multivariate_normal(
+ mean=mean_val, cov=cov_val, allow_singular=True
+ )
except ValueError as ve:
- logger.error(f"Scipy MvNormal error for {dist_name} with mean {mean_val} and cov {cov_val}: {ve}")
+ logger.error(
+ f"Scipy MvNormal error for {dist_name} with mean {mean_val} and cov {cov_val}: {ve}"
+ )
return None
return None
def _convert_mv_tdist():
# Distributions.jl: MvTDist(df, μ, Σ) where Σ is scale matrix (not covariance)
# scipy.stats.multivariate_t(loc=μ, shape=Σ, df=df) where shape is scale matrix
- df = _get_scalar_param("df") # degrees of freedom (from Distributions.jl MvTDist struct)
+ df = _get_scalar_param(
+ "df"
+ ) # degrees of freedom (from Distributions.jl MvTDist struct)
# Note: Distributions.jl MvTDist constructor is MvTDist(ν, μ, Σ), but fields might be df, loc, invscale, etc.
# Assuming direct field access for μ and Σ based on common parameter names.
# If Distributions.jl stores them as `loc` and `scale` for MvTDist, adjust "μ" and "Σ".
@@ -1155,21 +1261,32 @@ def distribution_from_julia(
if df is None:
df = _get_scalar_param("ν")
- mu_loc = _get_numpy_array_param("μ") # location vector
- sigma_shape = _get_numpy_array_param("Σ") # scale matrix
+ mu_loc = _get_numpy_array_param("μ") # location vector
+ sigma_shape = _get_numpy_array_param("Σ") # scale matrix
if df is not None and mu_loc is not None and sigma_shape is not None:
if df <= 0:
- logger.error(f"MvTDist 'df' (degrees of freedom) must be positive for {dist_name}, got {df}.")
+ logger.error(
+ f"MvTDist 'df' (degrees of freedom) must be positive for {dist_name}, got {df}."
+ )
return None
if mu_loc.ndim != 1:
- logger.error(f"MvTDist 'μ' (loc) must be 1D for {dist_name}, got {mu_loc.ndim}D.")
+ logger.error(
+ f"MvTDist 'μ' (loc) must be 1D for {dist_name}, got {mu_loc.ndim}D."
+ )
return None
if sigma_shape.ndim != 2:
- logger.error(f"MvTDist 'Σ' (shape) must be 2D for {dist_name}, got {sigma_shape.ndim}D.")
+ logger.error(
+ f"MvTDist 'Σ' (shape) must be 2D for {dist_name}, got {sigma_shape.ndim}D."
+ )
return None
- if mu_loc.shape[0] != sigma_shape.shape[0] or sigma_shape.shape[0] != sigma_shape.shape[1]:
- logger.error(f"MvTDist loc/shape dimension mismatch for {dist_name}: loc {mu_loc.shape}, shape {sigma_shape.shape}")
+ if (
+ mu_loc.shape[0] != sigma_shape.shape[0]
+ or sigma_shape.shape[0] != sigma_shape.shape[1]
+ ):
+ logger.error(
+ f"MvTDist loc/shape dimension mismatch for {dist_name}: loc {mu_loc.shape}, shape {sigma_shape.shape}"
+ )
return None
try:
# allow_singular might not be an option for multivariate_t, check scipy docs if issues arise
@@ -1180,17 +1297,21 @@ def distribution_from_julia(
return None
def _convert_dirichlet():
- alpha_conc = _get_numpy_array_param("alpha") # concentration parameters
+ alpha_conc = _get_numpy_array_param("alpha") # concentration parameters
# Distributions.jl uses 'α', let's try that first.
if alpha_conc is None:
alpha_conc = _get_numpy_array_param("α")
if alpha_conc is not None:
if alpha_conc.ndim != 1:
- logger.error(f"Dirichlet 'α' (concentration) must be 1D for {dist_name}, got {alpha_conc.ndim}D.")
+ logger.error(
+ f"Dirichlet 'α' (concentration) must be 1D for {dist_name}, got {alpha_conc.ndim}D."
+ )
return None
if not np.all(alpha_conc > 0):
- logger.error(f"Dirichlet 'α' elements must be positive for {dist_name}, got {alpha_conc}.")
+ logger.error(
+ f"Dirichlet 'α' elements must be positive for {dist_name}, got {alpha_conc}."
+ )
return None
return stats.dirichlet(alpha=alpha_conc)
return None
@@ -1199,7 +1320,9 @@ def distribution_from_julia(
p = _get_scalar_param("p")
if p is not None:
if not (0 <= p <= 1):
- logger.error(f"Bernoulli 'p' must be in [0,1] for {dist_name}, got {p}.")
+ logger.error(
+ f"Bernoulli 'p' must be in [0,1] for {dist_name}, got {p}."
+ )
return None
return stats.bernoulli(p=p)
return None
@@ -1210,13 +1333,19 @@ def distribution_from_julia(
p_probs = _get_numpy_array_param("p")
if p_probs is not None:
if p_probs.ndim != 1:
- logger.error(f"Categorical 'p' (probabilities) must be 1D for {dist_name}, got {p_probs.ndim}D.")
+ logger.error(
+ f"Categorical 'p' (probabilities) must be 1D for {dist_name}, got {p_probs.ndim}D."
+ )
return None
if not np.all(p_probs >= 0):
- logger.error(f"Categorical 'p' elements must be non-negative for {dist_name}, got {p_probs}.")
+ logger.error(
+ f"Categorical 'p' elements must be non-negative for {dist_name}, got {p_probs}."
+ )
return None
if not np.isclose(np.sum(p_probs), 1.0):
- logger.warning(f"Categorical 'p' elements for {dist_name} do not sum to 1 (sum={np.sum(p_probs)}). Normalizing for scipy.")
+ logger.warning(
+ f"Categorical 'p' elements for {dist_name} do not sum to 1 (sum={np.sum(p_probs)}). Normalizing for scipy."
+ )
p_probs = p_probs / np.sum(p_probs)
return stats.multinomial(n=1, p=p_probs)
return None
@@ -1225,21 +1354,43 @@ def distribution_from_julia(
n_trials = _get_int_param("n")
p_probs = _get_numpy_array_param("p")
if n_trials is not None and p_probs is not None:
- if n_trials < 0: # typically n_trials >= 1 for multinomial, but scipy allows n=0
- logger.error(f"Multinomial 'n' (trials) must be non-negative for {dist_name}, got {n_trials}.")
+ if (
+ n_trials < 0
+ ): # typically n_trials >= 1 for multinomial, but scipy allows n=0
+ logger.error(
+ f"Multinomial 'n' (trials) must be non-negative for {dist_name}, got {n_trials}."
+ )
return None
if p_probs.ndim != 1:
- logger.error(f"Multinomial 'p' (probabilities) must be 1D for {dist_name}, got {p_probs.ndim}D.")
+ logger.error(
+ f"Multinomial 'p' (probabilities) must be 1D for {dist_name}, got {p_probs.ndim}D."
+ )
return None
if not np.all(p_probs >= 0):
- logger.error(f"Multinomial 'p' elements must be non-negative for {dist_name}, got {p_probs}.")
+ logger.error(
+ f"Multinomial 'p' elements must be non-negative for {dist_name}, got {p_probs}."
+ )
return None
if not np.isclose(np.sum(p_probs), 1.0):
- logger.warning(f"Multinomial 'p' elements for {dist_name} do not sum to 1 (sum={np.sum(p_probs)}). Normalizing for scipy.")
+ logger.warning(
+ f"Multinomial 'p' elements for {dist_name} do not sum to 1 (sum={np.sum(p_probs)}). Normalizing for scipy."
+ )
p_probs = p_probs / np.sum(p_probs)
return stats.multinomial(n=n_trials, p=p_probs)
return None
+ def _convert_inversegamma1():
+ alpha_ig = _get_scalar_param("α")
+ theta_ig = _get_scalar_param("θ")
+ if alpha_ig is not None and theta_ig is not None:
+ if alpha_ig <= 0 or theta_ig <= 0:
+ logger.error(
+ f"InverseGamma1 'α' (shape) and 'θ' (scale) must be positive for {dist_name}, got α={alpha_ig}, θ={theta_ig}."
+ )
+ return None
+ return _inv_gamma1(a=alpha_ig, scale=theta_ig)
+ return None
+
# --- Conversion Map ---
conversion_map = {
"Normal": _convert_normal,
@@ -1251,8 +1402,8 @@ def distribution_from_julia(
"Poisson": _convert_poisson,
"Geometric": _convert_geometric,
"NegativeBinomial": _convert_negative_binomial,
- "ChiSquared": _convert_chi_squared, # Julia: ChiSquared, Scipy: chi2
- "Chisq": _convert_chi_squared, # Alias for ChiSquared
+ "ChiSquared": _convert_chi_squared, # Julia: ChiSquared, Scipy: chi2
+ "Chisq": _convert_chi_squared, # Alias for ChiSquared
"TDist": _convert_tdist,
"Laplace": _convert_laplace,
"Cauchy": _convert_cauchy,
@@ -1265,33 +1416,61 @@ def distribution_from_julia(
"Bernoulli": _convert_bernoulli,
"Categorical": _convert_categorical,
"Multinomial": _convert_multinomial,
- # Add InverseGamma from distributions.py if needed
- "InverseGamma": lambda: ( # Distributions.jl: InverseGamma(shape, scale) = InverseGamma(α, θ)
- alpha_ig := _get_scalar_param("α"), # shape
- theta_ig := _get_scalar_param("θ"), # scale
- (stats.invgamma(a=alpha_ig, scale=theta_ig)
- if alpha_ig is not None and theta_ig is not None and alpha_ig > 0 and theta_ig > 0
- else (logger.error(f"InverseGamma 'α' (shape) and 'θ' (scale) must be positive for {dist_name}, got α={alpha_ig}, θ={theta_ig}."), None)[1] if alpha_ig is not None and theta_ig is not None else None)
+ "InverseGamma": lambda: ( # Distributions.jl: InverseGamma(shape, scale) = InverseGamma(α, θ)
+ alpha_ig := _get_scalar_param("α"), # shape
+ theta_ig := _get_scalar_param("θ"), # scale
+ (
+ stats.invgamma(a=alpha_ig, scale=theta_ig)
+ if alpha_ig is not None
+ and theta_ig is not None
+ and alpha_ig > 0
+ and theta_ig > 0
+ else (
+ (
+ logger.error(
+ f"InverseGamma 'α' (shape) and 'θ' (scale) must be positive for {dist_name}, got α={alpha_ig}, θ={theta_ig}."
+ ),
+ None,
+ )[1]
+ if alpha_ig is not None and theta_ig is not None
+ else None
+ )
+ ),
),
- # Add FDist from distributions.py if needed
- "FDist": lambda: (
+ "InverseGamma1": _convert_inversegamma1,
+ "FDist": lambda: (
dfn := _get_scalar_param("ν1"),
dfd := _get_scalar_param("ν2"),
- (stats.f(dfn=dfn, dfd=dfd)
- if dfn is not None and dfd is not None and dfn > 0 and dfd > 0
- else (logger.error(f"FDist 'ν1' and 'ν2' must be positive for {dist_name}, got ν1={dfn}, ν2={dfd}."), None)[1] if dfn is not None and dfd is not None else None)
- )
+ (
+ stats.f(dfn=dfn, dfd=dfd)
+ if dfn is not None and dfd is not None and dfn > 0 and dfd > 0
+ else (
+ (
+ logger.error(
+ f"FDist 'ν1' and 'ν2' must be positive for {dist_name}, got ν1={dfn}, ν2={dfd}."
+ ),
+ None,
+ )[1]
+ if dfn is not None and dfd is not None
+ else None
+ )
+ ),
+ ),
}
if dist_name in conversion_map:
return conversion_map[dist_name]()
else:
- logger.warning(f"Unsupported Julia distribution type for scipy conversion: {dist_name}")
+ logger.warning(
+ f"Unsupported Julia distribution type for scipy conversion: {dist_name}"
+ )
return None
except Exception as e:
- # Use dist_name_for_log which might be "UnknownDistribution" or the actual name if obtained
- logger.error(f"Error converting Julia distribution ({dist_name_for_log}) to scipy: {e}", exc_info=True)
+ logger.error(
+ f"Error converting Julia distribution ({dist_name_for_log}) to scipy: {e}",
+ exc_info=True,
+ )
return None
@@ -1374,9 +1553,9 @@ class Work:
qr_jacobian: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
model_has_trend: List[bool] = field(default_factory=lambda: [False])
histval: np.ndarray = field(default_factory=lambda: np.empty((0, 0), dtype=object))
- homotopy_setup: List[
- Tuple[str, "SymbolType", int, float, Union[float, None]]
- ] = field(default_factory=list)
+ homotopy_setup: List[Tuple[str, "SymbolType", int, float, Union[float, None]]] = (
+ field(default_factory=list)
+ )
initval_endogenous: np.ndarray = field(
default_factory=lambda: np.empty((0, 0), dtype=object)
)
diff --git a/src/dynare/indices.py b/src/dynare/indices.py
index a4862b5..1f57359 100644
--- a/src/dynare/indices.py
+++ b/src/dynare/indices.py
@@ -1,6 +1,7 @@
from dataclasses import dataclass, field
from typing import List, Dict, Tuple, Any
+
@dataclass
class PyIndices:
current: List[int] = field(default_factory=list)
@@ -30,33 +31,49 @@ class PyIndices:
return cls()
d_cols_data = {}
- jl_D_columns = getattr(julia_indices, 'D_columns', None)
+ jl_D_columns = getattr(julia_indices, "D_columns", None)
if jl_D_columns is not None:
- if hasattr(jl_D_columns, 'D') and hasattr(jl_D_columns, 'jacobian'):
- d_cols_data = {
+ if hasattr(jl_D_columns, "D") and hasattr(jl_D_columns, "jacobian"):
+ d_cols_data = {
"D": list(jl_D_columns.D) if jl_D_columns.D is not None else [],
- "jacobian": list(jl_D_columns.jacobian) if jl_D_columns.jacobian is not None else [],
+ "jacobian": (
+ list(jl_D_columns.jacobian)
+ if jl_D_columns.jacobian is not None
+ else []
+ ),
}
- elif isinstance(jl_D_columns, tuple) and len(jl_D_columns) == 2: # Handles pre-converted tuple
- d_cols_data = {
+ elif (
+ isinstance(jl_D_columns, tuple) and len(jl_D_columns) == 2
+ ): # Handles pre-converted tuple
+ d_cols_data = {
"D": list(jl_D_columns[0]) if jl_D_columns[0] is not None else [],
- "jacobian": list(jl_D_columns[1]) if jl_D_columns[1] is not None else [],
+ "jacobian": (
+ list(jl_D_columns[1]) if jl_D_columns[1] is not None else []
+ ),
}
e_cols_data = {}
- jl_E_columns = getattr(julia_indices, 'E_columns', None)
+ jl_E_columns = getattr(julia_indices, "E_columns", None)
if jl_E_columns is not None:
- if hasattr(jl_E_columns, 'E') and hasattr(jl_E_columns, 'jacobian'):
+ if hasattr(jl_E_columns, "E") and hasattr(jl_E_columns, "jacobian"):
e_cols_data = {
"E": list(jl_E_columns.E) if jl_E_columns.E is not None else [],
- "jacobian": list(jl_E_columns.jacobian) if jl_E_columns.jacobian is not None else [],
+ "jacobian": (
+ list(jl_E_columns.jacobian)
+ if jl_E_columns.jacobian is not None
+ else []
+ ),
}
- elif isinstance(jl_E_columns, tuple) and len(jl_E_columns) == 2: # Handles pre-converted tuple
+ elif (
+ isinstance(jl_E_columns, tuple) and len(jl_E_columns) == 2
+ ): # Handles pre-converted tuple
e_cols_data = {
"E": list(jl_E_columns[0]) if jl_E_columns[0] is not None else [],
- "jacobian": list(jl_E_columns[1]) if jl_E_columns[1] is not None else [],
+ "jacobian": (
+ list(jl_E_columns[1]) if jl_E_columns[1] is not None else []
+ ),
}
-
+
# Helper to safely convert attributes that are expected to be lists
def _to_list(attr_name: str) -> List[Any]:
attr_val = getattr(julia_indices, attr_name, [])
@@ -64,30 +81,30 @@ class PyIndices:
return []
try:
return list(attr_val)
- except TypeError: # If it's not iterable (e.g. an int like n_endogenous)
+ except TypeError: # If it's not iterable (e.g. an int like n_endogenous)
return []
return cls(
- current=_to_list('current'),
- forward=_to_list('forward'),
- purely_forward=_to_list('purely_forward'),
- backward=_to_list('backward'),
- both=_to_list('both'),
- non_backward=_to_list('non_backward'),
- static=_to_list('static'),
- dynamic=_to_list('dynamic'),
- dynamic_current=_to_list('dynamic_current'),
- current_in_dynamic=_to_list('current_in_dynamic'),
- forward_in_dynamic=_to_list('forward_in_dynamic'),
- backward_in_dynamic=_to_list('backward_in_dynamic'),
- current_in_dynamic_jacobian=_to_list('current_in_dynamic_jacobian'),
- current_in_static_jacobian=_to_list('current_in_static_jacobian'),
- exogenous=_to_list('exogenous'),
- n_endogenous=getattr(julia_indices, 'n_endogenous', 0),
+ current=_to_list("current"),
+ forward=_to_list("forward"),
+ purely_forward=_to_list("purely_forward"),
+ backward=_to_list("backward"),
+ both=_to_list("both"),
+ non_backward=_to_list("non_backward"),
+ static=_to_list("static"),
+ dynamic=_to_list("dynamic"),
+ dynamic_current=_to_list("dynamic_current"),
+ current_in_dynamic=_to_list("current_in_dynamic"),
+ forward_in_dynamic=_to_list("forward_in_dynamic"),
+ backward_in_dynamic=_to_list("backward_in_dynamic"),
+ current_in_dynamic_jacobian=_to_list("current_in_dynamic_jacobian"),
+ current_in_static_jacobian=_to_list("current_in_static_jacobian"),
+ exogenous=_to_list("exogenous"),
+ n_endogenous=getattr(julia_indices, "n_endogenous", 0),
D_columns=d_cols_data,
E_columns=e_cols_data,
- UD_columns=_to_list('UD_columns'),
- UE_columns=_to_list('UE_columns'),
+ UD_columns=_to_list("UD_columns"),
+ UE_columns=_to_list("UE_columns"),
)
def get_D_columns_as_tuple(self) -> Tuple[List[int], List[int]]:
--
GitLab
From 210de44b82e8d303d62b1925c937c1b484ba6f15 Mon Sep 17 00:00:00 2001
From: Daniel Sali <daniel.sali@alphacruncher.com>
Date: Tue, 3 Jun 2025 09:30:34 +0200
Subject: [PATCH 3/6] Add docstring to dynare() function with example
---
examples/fs2000.mod | 168 +++++++++++++++++++++++++++++++++++++++++++
src/dynare/dynare.py | 26 +++++++
2 files changed, 194 insertions(+)
create mode 100644 examples/fs2000.mod
diff --git a/examples/fs2000.mod b/examples/fs2000.mod
new file mode 100644
index 0000000..13b95f0
--- /dev/null
+++ b/examples/fs2000.mod
@@ -0,0 +1,168 @@
+/*
+ * This file replicates the estimation of the cash in advance model (termed M1
+ * in the paper) described in Frank Schorfheide (2000): "Loss function-based
+ * evaluation of DSGE models", Journal of Applied Econometrics, 15(6), 645-670.
+ *
+ * The data are taken from the replication package at
+ * http://dx.doi.org/10.15456/jae.2022314.0708799949
+ *
+ * The prior distribution follows the one originally specified in Schorfheide's
+ * paper. Note that the elicited beta prior for rho in the paper
+ * implies an asymptote and corresponding prior mode at 0. It is generally
+ * recommended to avoid this extreme type of prior.
+ *
+ * Because the data are already logged and we use the loglinear option to conduct
+ * a full log-linearization, we need to use the logdata option.
+ *
+ * The equations are taken from J. Nason and T. Cogley (1994): "Testing the
+ * implications of long-run neutrality for monetary business cycle models",
+ * Journal of Applied Econometrics, 9, S37-S70, NC in the following.
+ * Note that there is an initial minus sign missing in equation (A1), p. S63.
+ *
+ * This implementation was originally written by Michel Juillard. Please note that the
+ * following copyright notice only applies to this Dynare implementation of the
+ * model.
+ */
+
+/*
+ * Copyright © 2004-2023 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+var m ${m}$ (long_name='money growth')
+ P ${P}$ (long_name='Price level')
+ c ${c}$ (long_name='consumption')
+ e ${e}$ (long_name='capital stock')
+ W ${W}$ (long_name='Wage rate')
+ R ${R}$ (long_name='interest rate')
+ k ${k}$ (long_name='capital stock')
+ d ${d}$ (long_name='dividends')
+ n ${n}$ (long_name='labor')
+ l ${l}$ (long_name='loans')
+ gy_obs ${\Delta \ln GDP}$ (long_name='detrended capital stock')
+ gp_obs ${\Delta \ln P}$ (long_name='detrended capital stock')
+ y ${y}$ (long_name='detrended output')
+ dA ${\Delta A}$ (long_name='TFP growth')
+ ;
+varexo e_a ${\epsilon_A}$ (long_name='TFP shock')
+ e_m ${\epsilon_M}$ (long_name='Money growth shock')
+ ;
+
+parameters alp ${\alpha}$ (long_name='capital share')
+ bet ${\beta}$ (long_name='discount factor')
+ gam ${\gamma}$ (long_name='long-run TFP growth')
+ logmst ${\log(m^*)}$ (long_name='long-run money growth')
+ rho ${\rho}$ (long_name='autocorrelation money growth')
+ phi ${\phi}$ (long_name='labor weight in consumption')
+ del ${\delta}$ (long_name='depreciation rate')
+ ;
+
+% roughly picked values to allow simulating the model before estimation
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+logmst = log(1.011);
+rho = 0.7;
+phi = 0.787;
+del = 0.02;
+
+model;
+[name='NC before eq. (1), TFP growth equation']
+dA = exp(gam+e_a);
+[name='NC eq. (2), money growth rate']
+log(m) = (1-rho)*logmst + rho*log(m(-1))+e_m;
+[name='NC eq. (A1), Euler equation']
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+[name='NC below eq. (A1), firm borrowing constraint']
+W = l/n;
+[name='NC eq. (A2), intratemporal labour market condition']
+-(phi/(1-phi))*(c*P/(1-n))+l/n = 0;
+[name='NC below eq. (A2), credit market clearing']
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+[name='NC eq. (A3), credit market optimality']
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+[name='NC eq. (18), aggregate resource constraint']
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+[name='NC eq. (19), money market condition']
+P*c = m;
+[name='NC eq. (20), credit market equilibrium condition']
+m-1+d = l;
+[name='Definition TFP shock']
+e = exp(e_a);
+[name='Implied by NC eq. (18), production function']
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+[name='Observation equation GDP growth']
+gy_obs = dA*y/y(-1);
+[name='Observation equation price level']
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = exp(logmst);
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/m )^(-1);
+ nust = phi*m^2/( (1-alp)*(1-phi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = phi*m*n/( (1-phi)*(1-n) );
+ c = m/P;
+ d = l - m + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = m/bet;
+ W = l/n;
+ ist = y-c;
+ q = 1 - d;
+
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+
+steady;
+check;
+
+% Table 1 of Schorfheide (2000)
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+logmst, normal_pdf, 0.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+phi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+estimation(order=1, datafile=fs2000_data, loglinear,logdata, mode_compute=4, mh_replic=20000, nodiagnostic, mh_nblocks=2, mh_jscale=0.8, mode_check);
+
+%uncomment the following lines to generate LaTeX-code of the model equations
+%write_latex_original_model(write_equation_tags);
+%collect_latex_files;
diff --git a/src/dynare/dynare.py b/src/dynare/dynare.py
index a097bab..91665c0 100644
--- a/src/dynare/dynare.py
+++ b/src/dynare/dynare.py
@@ -9,6 +9,32 @@ logger = logging.getLogger("dynare.dynare")
def dynare(model: str | Path, *dynare_options: str) -> Context:
+ """
+ Run a Dynare model using Julia and return a Python context object.
+ This function launches Julia, loads the Dynare package, and executes the specified
+ .mod file with any additional Dynare command-line options. It converts the resulting
+ Dynare Context into a Python‐friendly Context object containing model results,
+ workspaces, and simulation data.
+ Args:
+ model (str | pathlib.Path):
+ Path to the Dynare .mod file to be executed.
+ *dynare_options (str):
+ Additional command-line options to pass to Dynare. Each option should be a
+ string, e.g. "-DB=\"A string with space\"" or "-DA=true".
+ Returns:
+ Context:
+ A Context instance that wraps the Julia Dynare.Context, with results,
+ symbol tables, work data, and any simulations converted to Python objects.
+ Raises:
+ DynareError:
+ If Julia fails to produce a Dynare context or any unexpected error occurs
+ during execution.
+ Example:
+ ctx = dynare(
+ Path(__file__).parent / "examples" / "fs2000.mod",
+ "-DB=\"A string with space\"",
+ "-DA=true"
+ """
try:
if isinstance(model, str):
model = Path(model)
--
GitLab
From 24fb181476f2a48ad8ee06569aa91fda2c30be06 Mon Sep 17 00:00:00 2001
From: Daniel Sali <daniel.sali@alphacruncher.com>
Date: Tue, 3 Jun 2025 12:49:31 +0200
Subject: [PATCH 4/6] Workspace shell implementation
---
src/dynare/juliapkg.json | 8 +-
src/dynare/workspace.py | 265 +++++++++++++++++++++++++++++----------
2 files changed, 201 insertions(+), 72 deletions(-)
diff --git a/src/dynare/juliapkg.json b/src/dynare/juliapkg.json
index 9556035..80a6f78 100644
--- a/src/dynare/juliapkg.json
+++ b/src/dynare/juliapkg.json
@@ -7,19 +7,19 @@
},
"DataFrames": {
"uuid": "a93c6f00-e57d-5684-b7b6-d8193f3e46c0",
- "version": "1.7.0"
+ "version": "=1.7.0"
},
"Dynare": {
"uuid": "5203de40-99df-439e-afbc-014de65cb9ef",
- "version": "0.9.18"
+ "version": "=0.9.18"
},
"LinearRationalExpectations": {
"uuid": "45f42fbc-210c-4ecc-9452-59ec793b9bfd",
- "version": "0.5.7"
+ "version": "=0.5.7"
},
"Tables": {
"uuid": "bd369af6-aec1-5ad0-b16a-f7cc5008161c",
- "version": "1.12.0"
+ "version": "=1.12.0"
}
}
}
\ No newline at end of file
diff --git a/src/dynare/workspace.py b/src/dynare/workspace.py
index 2ab1c49..70a7675 100644
--- a/src/dynare/workspace.py
+++ b/src/dynare/workspace.py
@@ -1,7 +1,11 @@
+import logging
+from juliacall import Main as jl
from dataclasses import dataclass
-from typing import List, Union, NamedTuple, Any
+from typing import List, Union, Any, Optional
import numpy as np
-from .indices import PyIndices # Added import
+from .indices import PyIndices
+
+logger = logging.getLogger("dynare.workspace")
@dataclass
@@ -266,7 +270,7 @@ class PyVarianceWs:
Σ_ns_ns: np.ndarray
stationary_variables: List[bool]
nonstationary_ws: List[PyNonstationaryVarianceWs]
- lre_ws: PyLinearRationalExpectationsWs
+ lre_ws: "PyLinearRationalExpectationsWs"
lyapd_ws: PyLyapdWs
@classmethod
@@ -312,6 +316,42 @@ class PyCholeskyPivotedWs:
return cls(work=julia_obj.work.to_numpy(), piv=julia_obj.piv.to_numpy())
+@dataclass
+class PyGeneralizedSchurWs:
+ """Placeholder for Julia GeneralizedSchurWs."""
+
+ # This field is to accommodate direct instantiation patterns
+ _internal_data: Any = None
+
+ @classmethod
+ def from_julia(cls, julia_obj: Any) -> Optional["PyGeneralizedSchurWs"]:
+ """Converts a Julia GeneralizedSchurWs object to its Python equivalent."""
+ logger.warning(
+ "PyGeneralizedSchurWs.from_julia is a placeholder and not fully implemented. "
+ f"Received Julia object of type: {jl.nameof(jl.typeof(julia_obj)).__str__() if julia_obj is not None else 'None'}. "
+ "Returning empty instance."
+ )
+ return cls(_internal_data=julia_obj)
+
+
+@dataclass
+class PyCyclicReductionWs:
+ """Placeholder for Julia CyclicReductionWs."""
+
+ # Placeholder field to store the original Julia object
+ _julia_obj: Any = None
+
+ @classmethod
+ def from_julia(cls, julia_obj: Any) -> Optional["PyCyclicReductionWs"]:
+ """Converts a Julia CyclicReductionWs object to its Python equivalent."""
+ logger.warning(
+ "PyCyclicReductionWs.from_julia is a placeholder and not fully implemented. "
+ f"Received Julia object of type: {jl.nameof(jl.typeof(julia_obj)).__str__() if julia_obj is not None else 'None'}. "
+ "Returning empty instance."
+ )
+ return cls(_julia_obj=julia_obj)
+
+
@dataclass
class PyGsSolverWs:
tmp1: np.ndarray
@@ -320,7 +360,7 @@ class PyGsSolverWs:
g2: np.ndarray
luws1: PyLUWs
luws2: PyLUWs
- schurws: PyGeneralizedSchurWs
+ schurws: "PyGeneralizedSchurWs"
@classmethod
def from_julia(cls, d: Any, n1: int) -> "PyGsSolverWs":
@@ -334,64 +374,147 @@ class PyGsSolverWs:
Returns:
An instance of PyGsSolverWs.
"""
- n = d.shape[0]
+ logger.warning(
+ "PyGsSolverWs.from_julia is partially implemented. Using placeholder logic."
+ )
+
+ # Get matrix dimensions
+ n = (
+ d.shape[0] if hasattr(d, "shape") else n1 * 2
+ ) # Fallback if d doesn't have shape
n2 = n - n1
- tmp1 = np.empty((n1, n1)) # Create similar matrix for tmp1
- tmp2 = np.empty((n1, n1)) # Create similar matrix for tmp2
- g1 = np.empty((n1, n1)) # Create similar matrix for g1
- g2 = np.empty((n2, n1)) # Create similar matrix for g2
- # Conversion of luws1 and luws2 using LUWs equivalents
- luws1 = PyLUWs(np.empty((n1, n1))) # Example placeholder
- luws2 = PyLUWs(np.empty((n2, n2))) # Example placeholder
+ # Create empty matrices
+ tmp1 = np.empty((n1, n1))
+ tmp2 = np.empty((n1, n1))
+ g1 = np.empty((n1, n1))
+ g2 = np.empty((n2, n1))
+
+ # Create placeholder LUWs objects
+ luws1 = PyLUWs(ipiv=[0] * n1) # Empty ipiv list
+ luws2 = PyLUWs(ipiv=[0] * n2) # Empty ipiv list
+
+ # Create placeholder GeneralizedSchurWs
+ schurws = PyGeneralizedSchurWs.from_julia(d)
- # Conversion of schurws using GeneralizedSchurWs equivalents
- schurws = PyGeneralizedSchurWs(d.to_numpy()) # Convert d to numpy array
+ return cls(
+ tmp1=tmp1,
+ tmp2=tmp2,
+ g1=g1,
+ g2=g2,
+ luws1=luws1,
+ luws2=luws2,
+ schurws=schurws,
+ )
@dataclass
class PyLinearGsSolverWs:
solver_ws: "PyGsSolverWs"
ids: PyIndices
- d: List[float]
- e: List[float]
+ d: np.ndarray
+ e: np.ndarray
@classmethod
def from_julia(cls, julia_solver_ws: Any) -> "PyLinearGsSolverWs":
"""Converts a Julia LinearGsSolverWs object to its Python equivalent."""
- return cls(
- solver_ws=PyGsSolverWs.from_julia(
- julia_solver_ws.solver_ws
- ), # Assuming PyGsSolverWs is defined
- ids=PyIndices.from_julia(julia_solver_ws.ids),
- d=julia_solver_ws.d.to_numpy(), # Convert Julia Matrix to NumPy array
- e=julia_solver_ws.e.to_numpy(), # Convert Julia Matrix to NumPy array
+ logger.warning(
+ "PyLinearGsSolverWs.from_julia may use incomplete implementations. "
+ f"Received Julia object of type: {jl.nameof(jl.typeof(julia_solver_ws)).__str__() if julia_solver_ws is not None else 'None'}."
)
+ # Determine n1 from the shape of d
+ d_numpy = (
+ julia_solver_ws.d.to_numpy()
+ if hasattr(julia_solver_ws, "d")
+ else np.array([])
+ )
+ n1 = d_numpy.shape[0] if hasattr(d_numpy, "shape") else 1
+
+ try:
+ return cls(
+ solver_ws=PyGsSolverWs.from_julia(julia_solver_ws.solver_ws, n1),
+ ids=PyIndices.from_julia(julia_solver_ws.ids),
+ d=d_numpy,
+ e=(
+ julia_solver_ws.e.to_numpy()
+ if hasattr(julia_solver_ws, "e")
+ else np.array([])
+ ),
+ )
+ except Exception as e:
+ logger.error(f"Error in PyLinearGsSolverWs.from_julia: {e}")
+ # Return a minimal instance to avoid crashing
+ return cls(
+ solver_ws=PyGsSolverWs(
+ tmp1=np.array([]),
+ tmp2=np.array([]),
+ g1=np.array([]),
+ g2=np.array([]),
+ luws1=PyLUWs(ipiv=[]),
+ luws2=PyLUWs(ipiv=[]),
+ schurws=PyGeneralizedSchurWs(),
+ ),
+ ids=PyIndices(),
+ d=np.array([]),
+ e=np.array([]),
+ )
+
@dataclass
class PyLinearCyclicReductionWs:
solver_ws: PyCyclicReductionWs
ids: PyIndices
- a: List[float]
- b: List[float]
- c: List[float]
- x: List[float]
+ a: np.ndarray
+ b: np.ndarray
+ c: np.ndarray
+ x: np.ndarray
@classmethod
def from_julia(cls, julia_cyclic_ws: Any) -> "PyLinearCyclicReductionWs":
"""Converts a Julia LinearCyclicReductionWs object to its Python equivalent."""
- return cls(
- solver_ws=PyCyclicReductionWs.from_julia(
- julia_cyclic_ws.solver_ws
- ), # Assuming PyCyclicReductionWs is defined
- ids=PyIndices.from_julia(julia_cyclic_ws.ids),
- a=julia_cyclic_ws.a.to_numpy(),
- b=julia_cyclic_ws.b.to_numpy(),
- c=julia_cyclic_ws.c.to_numpy(),
- x=julia_cyclic_ws.x.to_numpy(),
+ logger.warning(
+ "PyLinearCyclicReductionWs.from_julia may use incomplete implementations. "
+ f"Received Julia object of type: {jl.nameof(jl.typeof(julia_cyclic_ws)).__str__() if julia_cyclic_ws is not None else 'None'}."
)
+ try:
+ return cls(
+ solver_ws=PyCyclicReductionWs.from_julia(julia_cyclic_ws.solver_ws),
+ ids=PyIndices.from_julia(julia_cyclic_ws.ids),
+ a=(
+ julia_cyclic_ws.a.to_numpy()
+ if hasattr(julia_cyclic_ws, "a")
+ else np.array([])
+ ),
+ b=(
+ julia_cyclic_ws.b.to_numpy()
+ if hasattr(julia_cyclic_ws, "b")
+ else np.array([])
+ ),
+ c=(
+ julia_cyclic_ws.c.to_numpy()
+ if hasattr(julia_cyclic_ws, "c")
+ else np.array([])
+ ),
+ x=(
+ julia_cyclic_ws.x.to_numpy()
+ if hasattr(julia_cyclic_ws, "x")
+ else np.array([])
+ ),
+ )
+ except Exception as e:
+ logger.error(f"Error in PyLinearCyclicReductionWs.from_julia: {e}")
+ # Return a minimal instance to avoid crashing
+ return cls(
+ solver_ws=PyCyclicReductionWs(),
+ ids=PyIndices(),
+ a=np.array([]),
+ b=np.array([]),
+ c=np.array([]),
+ x=np.array([]),
+ )
+
@dataclass
class PyLinearRationalExpectationsWs:
@@ -417,35 +540,41 @@ class PyLinearRationalExpectationsWs:
@classmethod
def from_julia(cls, julia_lre_ws: Any) -> "PyLinearRationalExpectationsWs":
"""Converts a Julia LinearRationalExpectationsWs object to its Python equivalent."""
- solver_ws = (
- PyLinearGsSolverWs.from_julia(julia_lre_ws.solver_ws)
- if type(julia_lre_ws.solver_ws) == "LinearGsSolverWs"
- else PyLinearCyclicReductionWs.from_julia(julia_lre_ws.solver_ws)
- )
-
- return cls(
- ids=PyIndices.from_julia(julia_lre_ws.ids),
- jacobian_static=julia_lre_ws.jacobian_static.to_numpy(),
- qr_ws=PyQRWs.from_julia(julia_lre_ws.qr_ws), # Assuming PyQRWs is defined
- ormqr_ws=PyQROrmWs.from_julia(
- julia_lre_ws.ormqr_ws
- ), # Assuming PyQROrmWs is defined
- solver_ws=solver_ws,
- A_s=julia_lre_ws.A_s.to_numpy(),
- C_s=julia_lre_ws.C_s.to_numpy(),
- Gy_forward=julia_lre_ws.Gy_forward.to_numpy(),
- Gy_dynamic=julia_lre_ws.Gy_dynamic.to_numpy(),
- temp=julia_lre_ws.temp.to_numpy(),
- AGplusB_backward=julia_lre_ws.AGplusB_backward.to_numpy(),
- jacobian_forward=julia_lre_ws.jacobian_forward.to_numpy(),
- jacobian_current=julia_lre_ws.jacobian_current.to_numpy(),
- b10=julia_lre_ws.b10.to_numpy(),
- b11=julia_lre_ws.b11.to_numpy(),
- AGplusB=julia_lre_ws.AGplusB.to_numpy(),
- linsolve_static_ws=PyLUWs.from_julia(
- julia_lre_ws.linsolve_static_ws
- ), # Assuming PyLUWs is defined
- AGplusB_linsolve_ws=PyLUWs.from_julia(
- julia_lre_ws.AGplusB_linsolve_ws
- ), # Assuming PyLUWs is defined
- )
+ try:
+ # Use juliacall to get the type name instead of Python's type()
+ julia_solver_ws_type_name = jl.nameof(
+ jl.typeof(julia_lre_ws.solver_ws)
+ ).__str__()
+
+ logger.info(f"Solver workspace type: {julia_solver_ws_type_name}")
+
+ # Determine which solver_ws implementation to use based on Julia type
+ if julia_solver_ws_type_name == "LinearGsSolverWs":
+ solver_ws = PyLinearGsSolverWs.from_julia(julia_lre_ws.solver_ws)
+ else:
+ solver_ws = PyLinearCyclicReductionWs.from_julia(julia_lre_ws.solver_ws)
+
+ return cls(
+ ids=PyIndices.from_julia(julia_lre_ws.ids),
+ jacobian_static=julia_lre_ws.jacobian_static.to_numpy(),
+ qr_ws=PyQRWs.from_julia(julia_lre_ws.qr_ws),
+ ormqr_ws=PyQROrmWs.from_julia(julia_lre_ws.ormqr_ws),
+ solver_ws=solver_ws,
+ A_s=julia_lre_ws.A_s.to_numpy(),
+ C_s=julia_lre_ws.C_s.to_numpy(),
+ Gy_forward=julia_lre_ws.Gy_forward.to_numpy(),
+ Gy_dynamic=julia_lre_ws.Gy_dynamic.to_numpy(),
+ temp=julia_lre_ws.temp.to_numpy(),
+ AGplusB_backward=julia_lre_ws.AGplusB_backward.to_numpy(),
+ jacobian_forward=julia_lre_ws.jacobian_forward.to_numpy(),
+ jacobian_current=julia_lre_ws.jacobian_current.to_numpy(),
+ b10=julia_lre_ws.b10.to_numpy(),
+ b11=julia_lre_ws.b11.to_numpy(),
+ AGplusB=julia_lre_ws.AGplusB.to_numpy(),
+ linsolve_static_ws=PyLUWs.from_julia(julia_lre_ws.linsolve_static_ws),
+ AGplusB_linsolve_ws=PyLUWs.from_julia(julia_lre_ws.AGplusB_linsolve_ws),
+ )
+ except Exception as e:
+ logger.error(f"Error in PyLinearRationalExpectationsWs.from_julia: {e}")
+ # Return None or a minimal instance to avoid crashing
+ return None
--
GitLab
From 0f0624b460ec875fa4457e66ff04e5a62df5a460 Mon Sep 17 00:00:00 2001
From: Daniel Sali <daniel.sali@alphacruncher.com>
Date: Wed, 4 Jun 2025 12:21:40 +0200
Subject: [PATCH 5/6] Lazy initialization of the Julia environment on dynare()
call
---
src/dynare/__init__.py | 38 ++++++++++++++++++++++++++++++++++----
1 file changed, 34 insertions(+), 4 deletions(-)
diff --git a/src/dynare/__init__.py b/src/dynare/__init__.py
index d57010b..96fc9cd 100644
--- a/src/dynare/__init__.py
+++ b/src/dynare/__init__.py
@@ -1,5 +1,6 @@
import logging
-import os # Added import
+import os
+import functools
from juliacall import Main as jl
from .dynare import dynare
@@ -42,6 +43,16 @@ def setup_logging(
def check_julia_and_dynare():
+ """
+ Check Julia installation and report on Dynare package status.
+
+ If you wish to use an existing conda environment with Julia and Dynare,
+ please set the environment variables appropriate for your setup as described at:
+ https://juliapy.github.io/PythonCall.jl/stable/pythoncall/#pythoncall-config
+
+ Julia configuration can be customized as described at:
+ https://juliapy.github.io/PythonCall.jl/stable/juliacall/#julia-config
+ """
# Check Julia
julia_path = jl.seval("Sys.BINDIR")
_ = jl.seval("Base.active_project()")
@@ -55,14 +66,33 @@ def check_julia_and_dynare():
_has_run = False
-def _run_once():
+def ensure_initialized():
+ """
+ Ensure Julia and Dynare packages are initialized.
+ This can be called explicitly by users who want to control when initialization happens.
+ """
global _has_run
if not _has_run:
check_julia_and_dynare()
jl.seval("using Pkg")
_has_run = True
+ logger.debug("Julia environment initialized")
+
+
+def lazy_init(func):
+ """
+ Decorator to ensure initialization before calling a function.
+ """
+
+ @functools.wraps(func)
+ def wrapper(*args, **kwargs):
+ ensure_initialized()
+ return func(*args, **kwargs)
+
+ return wrapper
-_run_once()
+# Make dynare function lazy-initialized
+dynare = lazy_init(dynare)
-__all__ = ["dynare", "setup_logging"]
+__all__ = ["dynare", "setup_logging", "ensure_initialized"]
--
GitLab
From fbe086dd34ebdcb8a441a187e1c67eb1848ae7ea Mon Sep 17 00:00:00 2001
From: Daniel Sali <daniel.sali@alphacruncher.com>
Date: Wed, 4 Jun 2025 12:57:55 +0200
Subject: [PATCH 6/6] Version bump to 0.2.0
---
pyproject.toml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/pyproject.toml b/pyproject.toml
index 49a6e4f..ac8ec28 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
[project]
name = "dynare-python"
-version = "0.1.0"
+version = "0.2.0"
description = "A Python wrapper for the Dynare.jl Julia package"
authors = [
{ name = "Daniel Sali", email = "daniel.sali@alphacruncher.com" }
--
GitLab