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