diff --git a/corrai/base/model.py b/corrai/base/model.py
index 551e8fa..9a82ca1 100644
--- a/corrai/base/model.py
+++ b/corrai/base/model.py
@@ -1,3 +1,4 @@
+import datetime as dt
 from typing import Union
 
 import pandas as pd
@@ -346,3 +347,39 @@ def simulate(
                 setattr(self, prop, val)
 
         return pd.Series({"res": self.prop_1 * self.prop_2 + self.prop_3})
+
+
+class Sine(PyModel):
+    def __init__(self):
+        super().__init__(is_dynamic=True)
+        self.omega = 2
+        self.amplitude = 5
+
+    def simulate(
+        self,
+        property_dict: dict[str, str | int | float] = None,
+        simulation_options: dict = None,
+        **simulation_kwargs,
+    ) -> pd.DataFrame | pd.Series:
+        self.set_property_values(property_dict)
+
+        start = simulation_options.get("start", "2009-01-01 00:00:00")
+        stop = simulation_options.get("stop", "2009-01-02 00:00:00")
+        output_freq = simulation_options.get("freq", "h")
+
+        index = pd.date_range(start, stop, freq=output_freq, tz="UTC")
+        cumsum_second = np.arange(
+            0, (index[-1] - index[0]).total_seconds() + 1, step=3600
+        )
+
+        return pd.DataFrame(
+            data=self.amplitude
+            * np.sin(
+                self.omega
+                * np.pi
+                / dt.timedelta(days=1).total_seconds()
+                * cumsum_second
+            ),
+            columns=["res"],
+            index=index,
+        )
diff --git a/corrai/base/utils.py b/corrai/base/utils.py
index 3f1ad80..a55113e 100644
--- a/corrai/base/utils.py
+++ b/corrai/base/utils.py
@@ -2,6 +2,7 @@
 import numpy as np
 import datetime as dt
 from collections.abc import Iterable
+from typing import Callable
 
 
 def _reshape_1d(sample):
@@ -258,3 +259,28 @@ def get_reversed_dict(dictionary, values=None):
         values = [values]
 
     return {val: key for key, val in dictionary.items() if val in values}
+
+
+def check_indicators_configs(
+    is_dynamic: bool,
+    indicators_configs: list[str]
+    | list[tuple[str, str | Callable] | tuple[str, str | Callable, pd.Series]]
+    | None,
+):
+    if is_dynamic:
+        if indicators_configs is None:
+            raise ValueError(
+                "Model is dynamic. At least one indicators and its aggregation "
+                "method must be provided"
+            )
+        if isinstance(indicators_configs[0], str):
+            raise ValueError(
+                "Invalid 'indicators_configs'. Model is dynamic"
+                "At least 'method' is required"
+            )
+    else:
+        if indicators_configs is not None and isinstance(indicators_configs[0], tuple):
+            raise ValueError(
+                "Invalid 'indicators_configs'. Model is static. "
+                "'indicators_configs' must be a list of string"
+            )
diff --git a/corrai/optimize.py b/corrai/optimize.py
index 93df99d..1b9a0ce 100644
--- a/corrai/optimize.py
+++ b/corrai/optimize.py
@@ -14,6 +14,7 @@
 
 from corrai.base.math import METHODS
 from corrai.base.model import Model
+from corrai.base.utils import check_indicators_configs
 from corrai.sampling import Sample
 from corrai.base.parameter import Parameter
 
@@ -209,35 +210,18 @@ def evaluate(
             np.array([[val[1] for val in parameter_value_pairs]]), [res]
         )
 
+        check_indicators_configs(self.model.is_dynamic, indicators_configs)
+
         if self.model.is_dynamic:
-            if indicators_configs is None:
-                raise ValueError(
-                    "Model is dynamic. At least one indicators and its aggregation "
-                    "method must be provided"
-                )
-            if isinstance(indicators_configs[0], str):
-                raise ValueError(
-                    "Invalid 'indicators_configs'. Model is dynamic"
-                    "At least 'method' is required"
-                )
             results = pd.Series()
             for config in indicators_configs:
                 col, func, *extra = config
                 series = res[col]
-
                 if isinstance(func, str):
                     func = METHODS[func]
-
                 results[col] = func(series, *extra)
             return pd.Series(results)
         else:
-            if indicators_configs is not None and isinstance(
-                indicators_configs[0], tuple
-            ):
-                raise ValueError(
-                    "Invalid 'indicators_configs'. Model is static. "
-                    "'indicators_configs' must be a list of string"
-                )
             return res[indicators_configs] if indicators_configs is not None else res
 
     def scipy_obj_function(self, x: np.ndarray, *args) -> float:
@@ -830,3 +814,48 @@ def plot_sample(
             quantile_band=quantile_band,
             type_graph=type_graph,
         )
+
+    @wraps(Sample.plot_pcp)
+    def plot_pcp(
+        self,
+        indicators_configs: list[str]
+        | list[tuple[str, str | Callable] | tuple[str, str | Callable, pd.Series]],
+        color_by: str | None = None,
+        title: str | None = "Parallel Coordinates — Samples",
+        html_file_path: str | None = None,
+    ) -> go.Figure:
+        return self.model_evaluator.sample.plot_pcp(
+            indicators_configs=indicators_configs,
+            color_by=color_by,
+            title=title,
+            html_file_path=html_file_path,
+        )
+
+    @wraps(Sample.plot_hist)
+    def plot_hist(
+        self,
+        indicator: str,
+        method: str = "mean",
+        unit: str = "",
+        agg_method_kwarg: dict = None,
+        reference_time_series: pd.Series = None,
+        bins: int = 30,
+        colors: str = "orange",
+        reference_value: int | float = None,
+        reference_label: str = "Reference",
+        show_rug: bool = False,
+        title: str = None,
+    ):
+        return self.model_evaluator.sample.plot_hist(
+            indicator=indicator,
+            method=method,
+            unit=unit,
+            agg_method_kwarg=agg_method_kwarg,
+            reference_time_series=reference_time_series,
+            bins=bins,
+            colors=colors,
+            reference_value=reference_value,
+            reference_label=reference_label,
+            show_rug=show_rug,
+            title=title,
+        )
diff --git a/corrai/sampling.py b/corrai/sampling.py
index c4f7a46..05de057 100644
--- a/corrai/sampling.py
+++ b/corrai/sampling.py
@@ -1,7 +1,7 @@
 from abc import ABC, abstractmethod
 from dataclasses import dataclass, field
 from functools import wraps
-from typing import Union
+from typing import Union, Callable
 
 import numpy as np
 import pandas as pd
@@ -15,6 +15,7 @@
 from SALib.sample import sobol as sobol_sampler
 from SALib.sample import fast_sampler, latin
 
+from corrai.base.utils import check_indicators_configs
 from corrai.base.parameter import Parameter
 from corrai.base.model import Model
 from corrai.base.math import aggregate_time_series
@@ -22,41 +23,34 @@
 
 
 def plot_pcp(
-    parameter_values: np.ndarray,
-    parameter_names: list[str],
+    parameter_values: pd.DataFrame,
     aggregated_results: pd.DataFrame,
-    *,
     bounds: list[tuple[float, float]] | None = None,
     color_by: str | None = None,
     title: str | None = "Parallel Coordinates — Samples",
     html_file_path: str | None = None,
 ) -> go.Figure:
     """
-    Creates a Parallel Coordinates Plot (PCP) for parameter samples and aggregated indicators.
-    Each vertical axis corresponds to a parameter or an aggregated indicator,
-    and each polyline represents one simulation.
+    Creates a Parallel Coordinates Plot (PCP) for parameter samples and aggregated
+    indicators. Each vertical axis corresponds to a parameter or an aggregated
+    indicator, and each polyline represents one simulation.
     """
 
-    if parameter_values.shape[0] != len(aggregated_results):
-        raise ValueError("Mismatch between number of samples and aggregated results.")
-    if len(parameter_names) != parameter_values.shape[1]:
+    if parameter_values.shape[0] != aggregated_results.shape[0]:
         raise ValueError(
-            "`parameter_names` length must match parameter_values.shape[1]."
+            "Shape mismatch between parameter_values and aggregated_results"
         )
 
-    df = pd.DataFrame(
-        parameter_values, columns=parameter_names, index=aggregated_results.index
-    )
-    df = pd.concat([df, aggregated_results], axis=1)
+    df = pd.concat([parameter_values, aggregated_results], axis=1)
 
     if color_by is None:
         if not aggregated_results.empty:
             color_by = aggregated_results.columns[0]
         else:
-            color_by = parameter_names[0]
+            color_by = parameter_values.columns[0]
 
     dimensions = []
-    for j, pname in enumerate(parameter_names):
+    for j, pname in enumerate(parameter_values.columns):
         dim = {"label": pname, "values": df[pname].to_numpy()}
         if bounds is not None:
             lb, ub = bounds[j]
@@ -712,6 +706,94 @@ def _legend_for(i: int) -> str:
         )
         return fig
 
+    def plot_pcp(
+        self,
+        indicators_configs: list[str]
+        | list[tuple[str, str | Callable] | tuple[str, str | Callable, pd.Series]],
+        color_by: str | None = None,
+        title: str | None = "Parallel Coordinates — Samples",
+        html_file_path: str | None = None,
+    ) -> go.Figure:
+        """
+        This method produces an interactive PCP visualization that allows comparison
+        of model parameters against aggregated indicators from simulation results.
+        It supports both dynamic and static models.
+
+        For dynamic models, the specified indicators are aggregated across time using
+        the provided functions (e.g., "mean", "sum", error metrics). For static models,
+        the indicators are taken directly from the stored results.
+
+        Parameters
+        ----------
+        indicators_configs : list of str or list of tuple
+            Configuration of indicators to include in the plot.
+
+            - For dynamic models, each element must be a tuple of the form:
+              ``(indicator_name, method)`` or
+              ``(indicator_name, method, reference_series)``.
+
+              Here:
+                * `indicator_name` : str
+                  Column name in the simulation results to aggregate.
+                * `method` : str or Callable
+                  Aggregation function or metric to apply.
+                * `reference_series` : pandas.Series, optional
+                  Reference time series required for error-based methods
+                  (e.g., mean absolute error).
+
+            - For static models, a simple list of indicator names (str) is sufficient.
+
+        color_by : str, optional
+            Name of a parameter or result column to use for coloring the PCP lines.
+            If None, all lines are plotted in the same color.
+
+        title : str, default="Parallel Coordinates — Samples"
+            Title of the plot.
+
+        html_file_path : str, optional
+            If provided, saves the interactive plot as an HTML file at the specified
+            path.
+
+        Returns
+        -------
+        plotly.graph_objects.Figure
+            The generated parallel coordinates figure. The figure can be displayed
+            interactively in a Jupyter notebook, web browser, or exported to HTML.
+
+        Raises
+        ------
+        ValueError
+            If the `indicators_configs` are incompatible with the model type
+            (dynamic vs static).
+
+        See Also
+        --------
+        get_aggregated_time_series :
+            For details on supported aggregation methods and how indicator values
+            are computed for dynamic models.
+        """
+
+        check_indicators_configs(self.is_dynamic, indicators_configs)
+
+        if self.is_dynamic:
+            results = pd.DataFrame()
+            for config in indicators_configs:
+                col, func, *extra = config
+                results[f"{func}_{col}"] = self.get_aggregated_time_series(
+                    col, func, reference_time_series=None if not extra else extra[0]
+                )
+        else:
+            results = self.get_static_results_as_df()[indicators_configs]
+
+        return plot_pcp(
+            parameter_values=self.values,
+            aggregated_results=results,
+            bounds=self.get_parameters_intervals().tolist(),
+            color_by=color_by,
+            title=title,
+            html_file_path=html_file_path,
+        )
+
 
 class Sampler(ABC):
     """
@@ -862,48 +944,17 @@ def get_sample_aggregated_time_series(
             indicator, method, agg_method_kwarg, reference_time_series, freq, prefix
         )
 
+    @wraps(Sample.plot_pcp)
     def plot_pcp(
         self,
-        indicator: str | None = None,
-        method: str | list[str] = "mean",
-        agg_method_kwarg: dict = None,
-        reference_time_series: pd.Series = None,
-        freq: str | pd.Timedelta | dt.timedelta = None,
-        prefix: str | None = None,
-        bounds: list[tuple[float, float]] | None = None,
+        indicators_configs: list[str]
+        | list[tuple[str, str | Callable] | tuple[str, str | Callable, pd.Series]],
         color_by: str | None = None,
         title: str | None = "Parallel Coordinates — Samples",
         html_file_path: str | None = None,
     ) -> go.Figure:
-        if indicator is None:
-            aggregated = pd.DataFrame(index=range(len(self.values)))
-        else:
-            methods = [method] if isinstance(method, str) else method
-            dfs = []
-            for m in methods:
-                this_prefix = prefix if prefix is not None else m
-                agg = aggregate_time_series(
-                    results=self.results,
-                    indicator=indicator,
-                    method=m,
-                    agg_method_kwarg=agg_method_kwarg,
-                    reference_time_series=reference_time_series,
-                    freq=freq,
-                    prefix=this_prefix,
-                )
-                if agg is not None and not agg.empty:
-                    dfs.append(agg)
-            aggregated = (
-                pd.concat(dfs, axis=1)
-                if dfs
-                else pd.DataFrame(index=range(len(self.values)))
-            )
-
-        return plot_pcp(
-            parameter_values=self.values,
-            parameter_names=[p.name for p in self.parameters],
-            aggregated_results=aggregated,
-            bounds=bounds,
+        return self.sample.plot_pcp(
+            indicators_configs=indicators_configs,
             color_by=color_by,
             title=title,
             html_file_path=html_file_path,
diff --git a/corrai/sensitivity.py b/corrai/sensitivity.py
index a46ad0b..55a1b1a 100644
--- a/corrai/sensitivity.py
+++ b/corrai/sensitivity.py
@@ -1,6 +1,7 @@
 import warnings
 from abc import ABC, abstractmethod
 from functools import wraps
+from typing import Callable
 
 import datetime as dt
 import numpy as np
@@ -398,27 +399,17 @@ def plot_sample(
             type_graph=type_graph,
         )
 
+    @wraps(Sample.plot_pcp)
     def plot_pcp(
         self,
-        indicator: str | None = None,
-        method: str = "mean",
-        agg_method_kwarg: dict = None,
-        reference_time_series: pd.Series = None,
-        freq: str | pd.Timedelta | dt.timedelta = None,
-        prefix: str | None = None,
-        bounds: list[tuple[float, float]] | None = None,
+        indicators_configs: list[str]
+        | list[tuple[str, str | Callable] | tuple[str, str | Callable, pd.Series]],
         color_by: str | None = None,
-        title: str | None = "Parallel Coordinates - Samples",
+        title: str | None = "Parallel Coordinates — Samples",
         html_file_path: str | None = None,
     ) -> go.Figure:
-        return self.sampler.plot_pcp(
-            indicator=indicator,
-            method=method,
-            agg_method_kwarg=agg_method_kwarg,
-            reference_time_series=reference_time_series,
-            freq=freq,
-            prefix=prefix,
-            bounds=bounds,
+        return self.sampler.sample.plot_pcp(
+            indicators_configs=indicators_configs,
             color_by=color_by,
             title=title,
             html_file_path=html_file_path,
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
index fceaf63..12ec109 100644
--- a/tests/test_sampling.py
+++ b/tests/test_sampling.py
@@ -4,15 +4,19 @@
 
 from corrai.base.parameter import Parameter
 from corrai.sampling import (
-    plot_pcp,
     LHSSampler,
     MorrisSampler,
     SobolSampler,
     Sample,
 )
 
-from corrai.base.math import aggregate_time_series
-from corrai.base.model import IshigamiDynamic, Ishigami, PymodelDynamic, PymodelStatic
+from corrai.base.model import (
+    IshigamiDynamic,
+    Ishigami,
+    PymodelDynamic,
+    PymodelStatic,
+    Sine,
+)
 
 import pytest
 
@@ -290,41 +294,34 @@ def test_plot_hist(self):
         assert len(hist.x) == len(sampler.results)
 
     def test_plot_pcp(self):
-        t = pd.date_range("2025-01-01 00:00:00", periods=2, freq="h")
-        df1 = pd.DataFrame({"res": [1.0, 2.0]}, index=t)
-        df2 = pd.DataFrame({"res": [3.0, 4.0]}, index=t)
-        df3 = pd.DataFrame({"res": [5.0, 6.0]}, index=t)
-        results = pd.Series([df1, df2, df3])
-
-        param_names = ["p1", "p2"]
-        param_values = np.array(
-            [
-                [1.1, 2.2],
-                [3.3, 4.4],
-                [5.5, 6.6],
-            ]
+        # Dynamic
+        lhs = LHSSampler(
+            parameters=[
+                Parameter("omega", (2, 8), model_property="omega"),
+                Parameter("amplitude", (1, 3), model_property="omega"),
+            ],
+            model=Sine(),
+            simulation_options={},
         )
 
-        agg_sum = aggregate_time_series(
-            results, indicator="res", method="sum", prefix="sum"
-        )
-        agg_mean = aggregate_time_series(
-            results, indicator="res", method="mean", prefix="mean"
-        )
-        aggregated = pd.concat([agg_sum, agg_mean], axis=1)
-
-        fig = plot_pcp(
-            parameter_values=param_values,
-            parameter_names=param_names,
-            aggregated_results=aggregated,
-            color_by="sum_res",
-            title="Parallel Coordinates — Samples",
+        lhs.add_sample(15)
+
+        lhs.sample.plot_pcp([("res", "mean"), ("res", "sum")])
+
+        # Static
+        lhs = LHSSampler(
+            parameters=[
+                Parameter("x1", (-3.14, 3.14), model_property="x"),
+                Parameter("x2", (-3.14, 3.14), model_property="x2"),
+                Parameter("x3", (-3.14, 3.14), model_property="x3"),
+            ],
+            model=Ishigami(),
+            simulation_options={},
         )
 
-        assert isinstance(fig, go.Figure)
-        assert len(fig.data) == 1
-        pc = fig.data[0]
-        np.testing.assert_allclose(pc.dimensions[0]["values"], [1.1, 3.3, 5.5])  # p1
+        lhs.add_sample(100)
+
+        lhs.sample.plot_pcp(["res"])
 
     def test_plot_pcp_in_sampler(self):
         sampler = LHSSampler(
@@ -335,7 +332,7 @@ def test_plot_pcp_in_sampler(self):
         sampler.add_sample(3, 42, simulate=True)
 
         fig = sampler.plot_pcp(
-            indicator="res",
+            indicators_configs=[("res", "mean")],
         )
         assert isinstance(fig, go.Figure)
         assert len(fig.data) == 1