diff --git a/examples/fs2000.mod b/examples/fs2000.mod new file mode 100644 index 0000000000000000000000000000000000000000..13b95f071485005113da68913c60c1ba04beb886 --- /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/pyproject.toml b/pyproject.toml index 49a6e4fd526295adae4d4ce2d51933ce6d9f1d1c..ac8ec28c42495f507a02bbf6720be855226e7b76 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" } diff --git a/src/dynare/__init__.py b/src/dynare/__init__.py index 327b6f513f03b14a0f97d2e54720b1bcaa2cd78b..96fc9cd772f1aea49e3bb4caf155c87e6e622408 100644 --- a/src/dynare/__init__.py +++ b/src/dynare/__init__.py @@ -1,20 +1,61 @@ -import subprocess import logging +import os +import functools 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 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") - julia_project = jl.seval("Base.active_project()") + _ = jl.seval("Base.active_project()") logger.info(f"Using julia at: {julia_path}") # Installed packages @@ -25,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"] \ No newline at end of file +__all__ = ["dynare", "setup_logging", "ensure_initialized"] diff --git a/src/dynare/distributions.py b/src/dynare/distributions.py deleted file mode 100644 index 3c41e07ebd7fed723b71cffc0373065147da4239..0000000000000000000000000000000000000000 --- 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 bb17fab2ff7bf4bacb4e9577ce6bca87daea49c1..91665c07d936d92698255a83354cf1c4065cd203 100644 --- a/src/dynare/dynare.py +++ b/src/dynare/dynare.py @@ -8,189 +8,222 @@ from .errors import DynareError logger = logging.getLogger("dynare.dynare") -def dynare(model: str | Path) -> Context: +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) 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("\\", "\\\\") - jl.seval( - "using DataFrames, Tables, PythonCall, Dynare, LinearRationalExpectations" - ) + resolved_path_str = str(model.resolve()) + escaped_path_for_julia_literal = resolved_path_str.replace( + "\\", "\\\\" + ).replace('"', '\\"') - # 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 + julia_options_str_parts = [] + for opt_str in dynare_options: + escaped_opt_for_julia_literal = opt_str.replace("\\", "\\\\").replace( + '"', '\\"' ) - end + julia_options_str_parts.append(f'"{escaped_opt_for_julia_literal}"') + + julia_options_cmd_part = " ".join(julia_options_str_parts) + + jl.seval( + """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 a6d5d7db8a1daff269b2ce419756543f88b9582c..aa2bd00aa3863074cc07e2d25c3e59989edc4fd4 100644 --- a/src/dynare/dynare_context.py +++ b/src/dynare/dynare_context.py @@ -3,10 +3,13 @@ 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 scipy.special import gamma as _gamma_func # for InverseGamma1 +from juliacall import Main as jl # Added import +from .indices import PyIndices logger = logging.getLogger("dynare.dynare_context") @@ -20,17 +23,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 +38,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 +104,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 +132,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 +328,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"), @@ -432,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"), @@ -444,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", [])), @@ -490,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"), @@ -567,7 +480,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 +525,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 +577,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 +630,10 @@ 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 +654,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 +664,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,9 +689,12 @@ 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") + 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() @@ -790,13 +709,35 @@ 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 +762,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 +822,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, @@ -915,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 @@ -943,7 +888,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 +923,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 +936,542 @@ 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. + """ + + 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 + + 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 + + # --- 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 + + 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, + "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, + "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 + ) + ), + ), + "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 + ) + ), + ), + } + + 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: + logger.error( + f"Error converting Julia distribution ({dist_name_for_log}) to scipy: {e}", + exc_info=True, ) + return None @dataclass @@ -1077,7 +1517,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), @@ -1112,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) ) @@ -1175,7 +1616,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 +1668,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 +1710,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 +1721,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 2dfb6a2c228603fa619ac013d1936f07c95faaab..a10ea8488d5b60473bccdd9f08ab4ec3247c5551 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 0000000000000000000000000000000000000000..1f57359846cbfd19e0c6efd0d8ddc87e0a436e46 --- /dev/null +++ b/src/dynare/indices.py @@ -0,0 +1,140 @@ +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 6464d8816e01b724be10da6dc4fcceffa6caf87f..80a6f7891c319393c41b3e01bf8f5993738ac66b 100644 --- a/src/dynare/juliapkg.json +++ b/src/dynare/juliapkg.json @@ -1,21 +1,25 @@ { "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" + "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 3ed2942828f26bb701c35ed67915122f1afb7ea1..70a767505e047e008159d4baadc8a064372d266a 100644 --- a/src/dynare/workspace.py +++ b/src/dynare/workspace.py @@ -1,61 +1,11 @@ -from dataclasses import dataclass, field -from typing import List, Any, Dict, Tuple, Union, Self, NamedTuple +import logging +from juliacall import Main as jl +from dataclasses import dataclass +from typing import List, Union, Any, Optional import numpy as np +from .indices import PyIndices - -@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), - ) +logger = logging.getLogger("dynare.workspace") @dataclass @@ -64,7 +14,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 +24,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 +35,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 +58,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 +77,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 +95,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 +120,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 +134,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 +148,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 +159,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 +172,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 +188,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 +209,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 +242,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(), @@ -306,11 +270,12 @@ class PyVarianceWs: Σ_ns_ns: np.ndarray stationary_variables: List[bool] nonstationary_ws: List[PyNonstationaryVarianceWs] - lre_ws: PyLinearRationalExpectationsWs + lre_ws: "PyLinearRationalExpectationsWs" 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 +300,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,10 +311,47 @@ 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()) +@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 @@ -357,66 +360,161 @@ class PyGsSolverWs: g2: np.ndarray luws1: PyLUWs luws2: PyLUWs - schurws: PyGeneralizedSchurWs + schurws: "PyGeneralizedSchurWs" @classmethod - def from_julia(cls, d, n1): - n = d.shape[0] + 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. + """ + 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)) - # Conversion of schurws using GeneralizedSchurWs equivalents - schurws = PyGeneralizedSchurWs(d.to_numpy()) # Convert d to numpy array + # 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) + + 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) -> "PyLinearGsSolverWs": - 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 + def from_julia(cls, julia_solver_ws: Any) -> "PyLinearGsSolverWs": + """Converts a Julia LinearGsSolverWs object to its Python equivalent.""" + 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) -> "PyLinearCyclicReductionWs": - 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(), + def from_julia(cls, julia_cyclic_ws: Any) -> "PyLinearCyclicReductionWs": + """Converts a Julia LinearCyclicReductionWs object to its Python equivalent.""" + 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: @@ -440,36 +538,43 @@ class PyLinearRationalExpectationsWs: AGplusB_linsolve_ws: "PyLUWs" @classmethod - def from_julia(cls, julia_lre_ws) -> Self: - 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 - ) + def from_julia(cls, julia_lre_ws: Any) -> "PyLinearRationalExpectationsWs": + """Converts a Julia LinearRationalExpectationsWs object to its Python equivalent.""" + 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