diff --git a/docs/source/nn.functional.rst b/docs/source/nn.functional.rst
index 5d6e90d4..b30e8206 100644
--- a/docs/source/nn.functional.rst
+++ b/docs/source/nn.functional.rst
@@ -110,6 +110,7 @@ Other Functions
     d2
     ncdf
     npdf
+    pl
     realized_variance
     realized_volatility
     svi_variance
diff --git a/pfhedge/nn/functional.py b/pfhedge/nn/functional.py
index ae9e5ad1..9ffcc542 100644
--- a/pfhedge/nn/functional.py
+++ b/pfhedge/nn/functional.py
@@ -1,5 +1,6 @@
 from math import ceil
 from math import pi as kPI
+from typing import List
 from typing import Optional
 from typing import Tuple
 from typing import Union
@@ -374,7 +375,6 @@ def realized_variance(input: Tensor, dt: TensorOrScalar) -> Tensor:
     Realized variance :math:`\sigma^2` of the stock price :math:`S` is defined as:
 
     .. math::
-
         \sigma^2 = \frac{1}{T - 1} \sum_{i = 1}^{T - 1}
         \frac{1}{dt} \log(S_{i + 1} / S_i)^2
 
@@ -419,28 +419,36 @@ def realized_volatility(input: Tensor, dt: Union[Tensor, float]) -> Tensor:
     return realized_variance(input, dt=dt).sqrt()
 
 
-def terminal_value(
+def pl(
     spot: Tensor,
     unit: Tensor,
-    cost: float = 0.0,
+    cost: Optional[List[float]] = None,
     payoff: Optional[Tensor] = None,
     deduct_first_cost: bool = True,
+    deduct_final_cost: bool = False,
 ) -> Tensor:
-    r"""Returns the terminal portfolio value.
+    r"""Returns the final profit and loss of hedging.
 
-    The terminal value of a hedger's portfolio is given by
+    For
+    hedging instruments indexed by :math:`h = 1, \dots, H` and
+    time steps :math:`i = 1, \dots, T`,
+    the final profit and loss is given by
 
     .. math::
-
         \text{PL}(Z, \delta, S) =
-        - Z
-        + \sum_{i = 0}^{T - 2} \delta_{i - 1} (S_{i} - S_{i - 1})
-        - c \sum_{i = 0}^{T - 1} |\delta_{i} - \delta_{i - 1}| S_{i}
+            - Z
+            + \sum_{h = 1}^{H} \sum_{t = 1}^{T} \left[
+                    \delta^{(h)}_{t - 1} (S^{(h)}_{t} - S^{(h)}_{t - 1})
+                    - c^{(h)} |\delta^{(h)}_{t} - \delta^{(h)}_{t - 1}| S^{(h)}_{t}
+                \right] ,
 
-    where :math:`Z` is the payoff of the derivative, :math:`T` is the number of
-    time steps, :math:`S` is the spot price, :math:`\delta` is the signed number
-    of shares held at each time step.
-    We define :math:`\delta_0 = 0` for notational convenience.
+    where
+    :math:`Z` is the payoff of the derivative.
+    For each hedging instrument,
+    :math:`\{S^{(h)}_t ; t = 1, \dots, T\}` is the spot price,
+    :math:`\{\delta^{(h)}_t ; t = 1, \dots, T\}` is the number of shares
+    held at each time step.
+    We define :math:`\delta^{(h)}_0 = 0` for notational convenience.
 
     A hedger sells the derivative to its customer and
     obliges to settle the payoff at maturity.
@@ -458,8 +466,8 @@ def terminal_value(
         spot (torch.Tensor): The spot price of the underlying asset :math:`S`.
         unit (torch.Tensor): The signed number of shares of the underlying asset
             :math:`\delta`.
-        cost (float, default=0.0): The proportional transaction cost rate of
-            the underlying asset :math:`c`.
+        cost (list[float], default=None): The proportional transaction cost rate of
+            the underlying assets.
         payoff (torch.Tensor, optional): The payoff of the derivative :math:`Z`.
         deduct_first_cost (bool, default=True): Whether to deduct the transaction
             cost of the stock at the first time step.
@@ -467,31 +475,57 @@ def terminal_value(
             equation of the terminal value.
 
     Shape:
-        - spot: :math:`(N, *, T)` where
-          :math:`T` is the number of time steps and
-          :math:`*` means any number of additional dimensions.
-        - unit: :math:`(N, *, T)`
-        - payoff: :math:`(N, *)`
-        - output: :math:`(N, *)`.
+        - spot: :math:`(N, H, T)` where
+          :math:`N` is the number of paths,
+          :math:`H` is the number of hedging instruments, and
+          :math:`T` is the number of time steps.
+        - unit: :math:`(N, H, T)`
+        - payoff: :math:`(N)`
+        - output: :math:`(N)`.
 
     Returns:
         torch.Tensor
     """
+    # TODO(simaki): Support deduct_final_cost=True
+    assert not deduct_final_cost, "not supported"
+
     if spot.size() != unit.size():
         raise RuntimeError(f"unmatched sizes: spot {spot.size()}, unit {unit.size()}")
-    if payoff is not None and spot.size()[:-1] != payoff.size():
-        raise RuntimeError(
-            f"unmatched sizes: spot {spot.size()}, payoff {payoff.size()}"
-        )
+    if payoff is not None:
+        if payoff.dim() != 1 or spot.size(0) != payoff.size(0):
+            raise RuntimeError(
+                f"unmatched sizes: spot {spot.size()}, payoff {payoff.size()}"
+            )
+
+    output = unit[..., :-1].mul(spot.diff(dim=-1)).sum(dim=(-2, -1))
 
-    value = unit[..., :-1].mul(spot.diff(dim=-1)).sum(-1)
-    value += -cost * unit.diff(dim=-1).abs().mul(spot[..., 1:]).sum(-1)
     if payoff is not None:
-        value -= payoff
-    if deduct_first_cost:
-        value -= cost * unit[..., 0].abs() * spot[..., 0]
+        output -= payoff
 
-    return value
+    if cost is not None:
+        c = torch.tensor(cost).unsqueeze(0).unsqueeze(-1)
+        output -= (spot[..., 1:] * unit.diff(dim=-1).abs() * c).sum(dim=(-2, -1))
+        if deduct_first_cost:
+            output -= (spot[..., [0]] * unit[..., [0]].abs() * c).sum(dim=(-2, -1))
+
+    return output
+
+
+def terminal_value(
+    spot: Tensor,
+    unit: Tensor,
+    cost: Optional[List[float]] = None,
+    payoff: Optional[Tensor] = None,
+    deduct_first_cost: bool = True,
+) -> Tensor:
+    """Alias for :func:`pfhedge.nn.functional.pl`."""
+    return pl(
+        spot=spot,
+        unit=unit,
+        cost=cost,
+        payoff=payoff,
+        deduct_first_cost=deduct_first_cost,
+    )
 
 
 def ncdf(input: Tensor) -> Tensor:
@@ -727,7 +761,29 @@ def bs_european_price(
 ) -> Tensor:
     """Returns Black-Scholes price of a European option.
 
-    See :func:`pfhedge.nn.BSEuropeanOption.price` for details.
+    .. seealso::
+        - :func:`pfhedge.nn.BSEuropeanOption.price`
+
+    Args:
+        log_moneyness (torch.Tensor, optional): Log moneyness of the underlying asset.
+        time_to_maturity (torch.Tensor, optional): Time to expiry of the option.
+        volatility (torch.Tensor, optional): Volatility of the underlying asset.
+
+    Shape:
+        - log_moneyness: :math:`(N, *)` where
+          :math:`*` means any number of additional dimensions.
+        - time_to_maturity: :math:`(N, *)`
+        - volatility: :math:`(N, *)`
+        - output: :math:`(N, *)`
+
+    Returns:
+        torch.Tensor
+
+    Examples:
+        >>> from pfhedge.nn.functional import bs_european_price
+        ...
+        >>> bs_european_price(torch.tensor([-0.1, 0.0, 0.1]), 1.0, 0.2)
+        tensor([0.0375, 0.0797, 0.1467])
     """
     s, t, v = broadcast_all(log_moneyness, time_to_maturity, volatility)
 
@@ -746,7 +802,29 @@ def bs_european_delta(
 ) -> Tensor:
     """Returns Black-Scholes delta of a European option.
 
-    See :func:`pfhedge.nn.BSEuropeanOption.delta` for details.
+    .. seealso::
+        - :func:`pfhedge.nn.BSEuropeanOption.delta`
+
+    Args:
+        log_moneyness (torch.Tensor, optional): Log moneyness of the underlying asset.
+        time_to_maturity (torch.Tensor, optional): Time to expiry of the option.
+        volatility (torch.Tensor, optional): Volatility of the underlying asset.
+
+    Shape:
+        - log_moneyness: :math:`(N, *)` where
+          :math:`*` means any number of additional dimensions.
+        - time_to_maturity: :math:`(N, *)`
+        - volatility: :math:`(N, *)`
+        - output: :math:`(N, *)`
+
+    Returns:
+        torch.Tensor
+
+    Examples:
+        >>> from pfhedge.nn.functional import bs_european_delta
+        ...
+        >>> bs_european_delta(torch.tensor([-0.1, 0.0, 0.1]), 1.0, 0.2)
+        tensor([0.3446, 0.5398, 0.7257])
     """
     s, t, v = broadcast_all(log_moneyness, time_to_maturity, volatility)
 
@@ -764,7 +842,29 @@ def bs_european_gamma(
 ) -> Tensor:
     """Returns Black-Scholes gamma of a European option.
 
-    See :func:`pfhedge.nn.BSEuropeanOption.gamma` for details.
+    .. seealso::
+        - :func:`pfhedge.nn.BSEuropeanOption.gamma`
+
+    Args:
+        log_moneyness (torch.Tensor, optional): Log moneyness of the underlying asset.
+        time_to_maturity (torch.Tensor, optional): Time to expiry of the option.
+        volatility (torch.Tensor, optional): Volatility of the underlying asset.
+
+    Shape:
+        - log_moneyness: :math:`(N, *)` where
+          :math:`*` means any number of additional dimensions.
+        - time_to_maturity: :math:`(N, *)`
+        - volatility: :math:`(N, *)`
+        - output: :math:`(N, *)`
+
+    Returns:
+        torch.Tensor
+
+    Examples:
+        >>> from pfhedge.nn.functional import bs_european_gamma
+        ...
+        >>> bs_european_gamma(torch.tensor([-0.1, 0.0, 0.1]), 1.0, 0.2)
+        tensor([2.0350, 1.9848, 1.5076])
     """
     s, t, v = broadcast_all(log_moneyness, time_to_maturity, volatility)
     spot = strike * s.exp()
diff --git a/pfhedge/nn/modules/bs/lookback.py b/pfhedge/nn/modules/bs/lookback.py
index 3c868a11..fef86331 100644
--- a/pfhedge/nn/modules/bs/lookback.py
+++ b/pfhedge/nn/modules/bs/lookback.py
@@ -268,7 +268,7 @@ def gamma(
             torch.Tensor
 
         Note:
-            args are not optional if it doesn't accept derivative in this initialization.
+            Arguments are not optional if it doesn't accept derivative in this initialization.
         """
         (
             log_moneyness,
@@ -319,7 +319,7 @@ def vega(
             torch.Tensor
 
         Note:
-            args are not optional if it doesn't accept derivative in this initialization.
+            - Arguments are not optional if it doesn't accept derivative in this initialization.
         """
         (
             log_moneyness,
@@ -366,14 +366,12 @@ def theta(
             - volatility: :math:`(N, *)`
             - output: :math:`(N, *)`
 
-        Note:
-            Risk-free rate is set to zero.
-
         Returns:
             torch.Tensor
 
         Note:
-            args are not optional if it doesn't accept derivative in this initialization.
+            - Risk-free rate is set to zero.
+            - Arguments are not optional if it doesn't accept derivative in this initialization.
         """
         (
             log_moneyness,
@@ -424,9 +422,9 @@ def implied_volatility(
             torch.Tensor
 
         Note:
-            args are not optional if it doesn't accept derivative in this initialization.
-            price seems optional in typing, but it isn't. It is set for the compatibility to the previous versions.
+            - Arguments are not optional if it doesn't accept derivative in this initialization.
         """
+        # price seems optional in typing, but it isn't. It is set for the compatibility to the previous versions.
         (log_moneyness, time_to_maturity) = acquire_params_from_derivative_0(
             self.derivative, log_moneyness, time_to_maturity
         )
diff --git a/pfhedge/nn/modules/hedger.py b/pfhedge/nn/modules/hedger.py
index edb050a3..804b206c 100644
--- a/pfhedge/nn/modules/hedger.py
+++ b/pfhedge/nn/modules/hedger.py
@@ -24,7 +24,7 @@
 from pfhedge.features._base import Feature
 from pfhedge.instruments.base import BaseInstrument
 from pfhedge.instruments.derivative.base import BaseDerivative
-from pfhedge.nn.functional import terminal_value
+from pfhedge.nn.functional import pl
 
 from .loss import EntropicRiskMeasure
 from .loss import HedgeLoss
@@ -62,7 +62,6 @@ class Hedger(Module):
           :math:`H` is the number of hedging instruments.
 
     Examples:
-
         A hedger that uses Black-Scholes' delta hedging strategy.
         See :class:`pfhedge.nn.BlackScholes` for details of the module.
 
@@ -70,7 +69,7 @@ class Hedger(Module):
         >>> from pfhedge.instruments import EuropeanOption
         >>> from pfhedge.nn import BlackScholes
         >>> from pfhedge.nn import Hedger
-        >>>
+        ...
         >>> derivative = EuropeanOption(BrownianStock(cost=1e-4))
         >>> model = BlackScholes(derivative)
         >>> hedger = Hedger(model, model.inputs())
@@ -111,7 +110,8 @@ class Hedger(Module):
         >>>
         >>> model = MultiLayerPerceptron()
         >>> hedger = Hedger(model, ["moneyness", "time_to_maturity", "volatility"])
-        >>> _ = hedger.compute_pnl(derivative, n_paths=1)  # Lazily materialize
+        >>> derivative.simulate(n_paths=1)
+        >>> _ = hedger.compute_pl(derivative)  # Lazily materialize
         >>> hedger
         Hedger(
           inputs=['moneyness', 'time_to_maturity', 'volatility']
@@ -188,7 +188,7 @@ def __init__(
         model: Module,
         inputs: List[Union[str, Feature]],
         criterion: HedgeLoss = EntropicRiskMeasure(),
-    ):
+    ) -> None:
         super().__init__()
 
         self.model = model
@@ -233,11 +233,10 @@ def get_input(self, derivative: BaseDerivative, time_step: Optional[int]) -> Ten
             torch.Tensor
 
         Examples:
-
             >>> from pfhedge.instruments import BrownianStock
             >>> from pfhedge.instruments import EuropeanOption
             >>> from pfhedge.nn import Naked
-            >>>
+            ...
             >>> derivative = EuropeanOption(BrownianStock())
             >>> derivative.simulate()
             >>> hedger = Hedger(Naked(), ["time_to_maturity", "volatility"])
@@ -247,19 +246,27 @@ def get_input(self, derivative: BaseDerivative, time_step: Optional[int]) -> Ten
         """
         return self.inputs.of(derivative=derivative).get(time_step)
 
+    def _get_hedge(
+        self, derivative: BaseDerivative, hedge: Optional[List[BaseInstrument]]
+    ) -> List[BaseInstrument]:
+        if hedge is None:
+            hedge = list(derivative.underliers())
+        return cast(List[BaseInstrument], hedge)
+
     def compute_hedge(
         self, derivative: BaseDerivative, hedge: Optional[List[BaseInstrument]] = None
     ) -> Tensor:
         """Compute the hedge ratio at each time step.
-        It assumes that the derivative is already simulated.
+
+        This method assumes that the derivative is already simulated.
 
         Args:
             derivative (BaseDerivative): The derivative to hedge.
-            hedge (BaseInstrument, optional): The hedging instrument.
-                If ``None`` (default), use ``derivative.underlier``.
+            hedge (list[BaseInstrument], optional): The hedging instruments.
+                If ``None`` (default), use ``derivative.underliers``.
 
         Shape:
-            - Output: :math:`(N, H, T)` where
+            - output: :math:`(N, H, T)` where
               :math:`N` is the number of paths,
               :math:`H` is the number of hedging instruments, and
               :math:`T` is the number of time steps.
@@ -268,11 +275,10 @@ def compute_hedge(
             torch.Tensor
 
         Examples:
-
             >>> from pfhedge.instruments import BrownianStock
             >>> from pfhedge.instruments import EuropeanOption
             >>> from pfhedge.nn import BlackScholes
-            >>>
+            ...
             >>> _ = torch.manual_seed(42)
             >>> derivative = EuropeanOption(BrownianStock(), maturity=5/250)
             >>> derivative.simulate(n_paths=2)
@@ -286,9 +292,7 @@ def compute_hedge(
                     [0.5056, 0.3785, 0.4609, 0.5239, 0.7281, 0.7281]])
         """
         inputs = self.inputs.of(derivative, self)
-        hedge = cast(
-            List[BaseInstrument], [derivative.ul()] if hedge is None else hedge
-        )
+        hedge = self._get_hedge(derivative, hedge)
 
         # Check that the spot prices of the hedges have the same sizes
         if not all(h.spot.size() == hedge[0].spot.size() for h in hedge):
@@ -318,19 +322,44 @@ def compute_hedge(
 
         return output
 
-    def compute_pnl(
-        self,
-        derivative: BaseDerivative,
-        hedge: Optional[List[BaseInstrument]] = None,
-        n_paths: int = 1000,
-        init_state: Optional[Tuple[TensorOrScalar, ...]] = None,
+    def compute_portfolio(
+        self, derivative: BaseDerivative, hedge: Optional[List[BaseInstrument]] = None
+    ) -> Tensor:
+        r"""Compute terminal value of the hedging portfolio.
+
+        See :func:`pfhedge.nn.functional.pl`, with :math:`Z` being substituted with 0,
+        for the expression of the terminal value of the hedging portfolio.
+
+        This method assumes that the derivative is already simulated.
+
+        Args:
+            derivative (BaseDerivative): The derivative to hedge.
+            hedge (BaseInstrument, optional): The hedging instrument.
+                If ``None`` (default), use ``derivative.underlier``.
+
+        Shape:
+            - output: :math:`(N)` where :math:`N` is the number of paths.
+
+        Returns:
+            torch.Tensor
+        """
+        hedge = self._get_hedge(derivative, hedge)
+
+        spot = torch.stack([h.spot for h in hedge], dim=1)
+        unit = self.compute_hedge(derivative, hedge=hedge)
+        cost = [h.cost for h in hedge]
+
+        return pl(spot=spot, unit=unit, cost=cost)
+
+    def compute_pl(
+        self, derivative: BaseDerivative, hedge: Optional[List[BaseInstrument]] = None
     ) -> Tensor:
         """Returns the terminal portfolio value after hedging a given derivative.
 
-        This method simulates the derivative, computes the hedge ratio, and
-        computes the terminal portfolio value.
+        This method assumes that the derivative is already simulated.
+
         See :func:`pfhedge.nn.functional.terminal_value` for the expression of the
-        terminal portyol value after hedging a derivative.
+        terminal portfolio value after hedging a derivative.
 
         Args:
             derivative (BaseDerivative): The derivative to hedge.
@@ -350,30 +379,37 @@ def compute_pnl(
             torch.Tensor
 
         Examples:
-
             >>> from pfhedge.instruments import BrownianStock
             >>> from pfhedge.instruments import EuropeanOption
             >>> from pfhedge.nn import BlackScholes
             >>> from pfhedge.nn import Hedger
-            >>>
+            ...
             >>> derivative = EuropeanOption(BrownianStock())
+            >>> derivative.simulate(n_paths=2)
             >>> model = BlackScholes(derivative)
             >>> hedger = Hedger(model, model.inputs())
-            >>> hedger.compute_pnl(derivative, n_paths=2)
+            >>> hedger.compute_pl(derivative)
             tensor([..., ...])
         """
-        derivative.simulate(n_paths=n_paths, init_state=init_state)
-        hedge = cast(
-            List[BaseInstrument], [derivative.ul()] if hedge is None else hedge
-        )
+        hedge = self._get_hedge(derivative, hedge)
 
+        spot = torch.stack([h.spot for h in hedge], dim=1)
         unit = self.compute_hedge(derivative, hedge=hedge)
+        cost = [h.cost for h in hedge]
 
-        output = -derivative.payoff()
-        for i, h in enumerate(hedge):
-            output += terminal_value(h.spot, unit=unit[:, i, :], cost=h.cost)
+        return pl(spot=spot, unit=unit, cost=cost, payoff=derivative.payoff())
 
-        return output
+    def compute_pnl(
+        self,
+        derivative: BaseDerivative,
+        hedge: Optional[List[BaseInstrument]] = None,
+        n_paths: int = 1000,
+        init_state: Optional[Tuple[TensorOrScalar, ...]] = None,
+    ) -> Tensor:
+        """(deprecated) Simulates derivative and computes profit loss by :meth:`compute_pl`."""
+        # TODO(simaki): Raise DeprecationWarning later
+        derivative.simulate(n_paths=n_paths, init_state=init_state)
+        return self.compute_pl(derivative=derivative, hedge=hedge)
 
     def compute_loss(
         self,
@@ -387,8 +423,8 @@ def compute_loss(
         """Returns the value of the criterion for the terminal portfolio value
         after hedging a given derivative.
 
-        This method basically computes ``self.criterion(pnl)``
-        where ``pnl`` is given by :meth:`compute_pnl`.
+        This method basically computes ``self.criterion(pl)``
+        where ``pl`` is given by :meth:`compute_pl`.
 
         Args:
             derivative (BaseDerivative): The derivative to hedge.
@@ -412,25 +448,40 @@ def compute_loss(
             torch.Tensor
 
         Examples:
-
             >>> from pfhedge.instruments import BrownianStock
             >>> from pfhedge.instruments import EuropeanOption
             >>> from pfhedge.nn import BlackScholes
             >>> from pfhedge.nn import Hedger
-            >>>
+            ...
             >>> derivative = EuropeanOption(BrownianStock())
             >>> model = BlackScholes(derivative)
             >>> hedger = Hedger(model, model.inputs())
             >>> hedger.compute_loss(derivative, n_paths=2)
             tensor(...)
+
+            One can use PyTorch built-in loss functions,
+            such as the mean squared loss :class:`torch.nn.MSELoss`, as criteria.
+            Then the criterion measures the loss between the hedging portfolio
+            (cf. :meth:`compute_portfolio`) as ``input`` and
+            the payoff of the derivative as ``target``.
+
+            >>> from torch.nn import MSELoss
+            ...
+            >>> _ = torch.manual_seed(42)
+            >>> derivative = EuropeanOption(BrownianStock())
+            >>> model = BlackScholes(derivative)
+            >>> hedger = Hedger(model, model.inputs(), criterion=MSELoss())
+            >>> hedger.compute_loss(derivative, n_paths=10)
+            tensor(...)
         """
         with torch.set_grad_enabled(enable_grad):
-            loss = lambda: self.criterion(
-                self.compute_pnl(
-                    derivative, hedge=hedge, n_paths=n_paths, init_state=init_state
-                )
-            )
-            mean_loss = ensemble_mean(loss, n_times=n_times)
+
+            def _get_loss():
+                derivative.simulate(n_paths=n_paths, init_state=init_state)
+                portfolio = self.compute_portfolio(derivative, hedge=hedge)
+                return self.criterion(portfolio, derivative.payoff())
+
+            mean_loss = ensemble_mean(_get_loss, n_times=n_times)
 
         return mean_loss
 
@@ -442,7 +493,8 @@ def _configure_optimizer(
         if not isinstance(optimizer, Optimizer):
             if has_lazy(self):
                 # Run a placeholder forward to initialize lazy parameters
-                _ = self.compute_pnl(derivative, n_paths=1)
+                derivative.simulate(n_paths=1)
+                _ = self.compute_pl(derivative)
             # If we use `if issubclass(optimizer, Optimizer)` here, mypy thinks that
             # optimizer is Optimizer rather than its subclass (e.g. Adam)
             # and complains that the required parameter default is missing.
@@ -467,8 +519,8 @@ def fit(
     ) -> Optional[List[float]]:
         """Fit the hedging model to hedge a given derivative.
 
-        The training is performed so that the hedger minimizes ``criterion(pnl)``
-        where ``pnl`` is given by :meth:`compute_pnl`.
+        The training is performed so that the hedger minimizes ``criterion(pl)``
+        where ``pl`` is given by :meth:`compute_pl`.
 
         It returns the training history, that is,
         validation loss after each simulation.
@@ -499,11 +551,10 @@ def fit(
             list[float]
 
         Examples:
-
             >>> from pfhedge.instruments import BrownianStock
             >>> from pfhedge.instruments import EuropeanOption
             >>> from pfhedge.nn import MultiLayerPerceptron
-            >>>
+            ...
             >>> derivative = EuropeanOption(BrownianStock())
             >>> model = MultiLayerPerceptron()
             >>> hedger = Hedger(model, ["moneyness", "time_to_maturity", "volatility"])
@@ -515,7 +566,7 @@ def fit(
             >>> from pfhedge.instruments import EuropeanOption
             >>> from pfhedge.nn import MultiLayerPerceptron
             >>> from torch.optim import SGD
-            >>>
+            ...
             >>> derivative = EuropeanOption(BrownianStock())
             >>> hedger = Hedger(MultiLayerPerceptron(), ["empty"])
             >>> # Run a placeholder forward to initialize lazy parameters
@@ -530,7 +581,7 @@ def fit(
             The optimizer will be initialized as ``Adadelta(hedger.parameters())``.
 
             >>> from torch.optim import Adadelta
-            >>>
+            ...
             >>> derivative = EuropeanOption(BrownianStock())
             >>> hedger = Hedger(MultiLayerPerceptron(), ["empty"])
             >>> _ = hedger.fit(
@@ -603,12 +654,11 @@ def price(
             torch.Tensor
 
         Examples:
-
             >>> from pfhedge.instruments import BrownianStock
             >>> from pfhedge.instruments import EuropeanOption
             >>> from pfhedge.nn import BlackScholes
             >>> from pfhedge.nn import Hedger
-            >>>
+            ...
             >>> derivative = EuropeanOption(BrownianStock())
             >>> model = BlackScholes(derivative)
             >>> hedger = Hedger(model, model.inputs())
@@ -616,12 +666,13 @@ def price(
             tensor(...)
         """
         with torch.set_grad_enabled(enable_grad):
-            # Negative because selling
-            pricer = lambda: -self.criterion.cash(
-                self.compute_pnl(
-                    derivative, hedge=hedge, n_paths=n_paths, init_state=init_state
-                )
-            )
-            mean_price = ensemble_mean(pricer, n_times=n_times)
+
+            def _get_price():
+                derivative.simulate(n_paths=n_paths, init_state=init_state)
+                portfolio = self.compute_portfolio(derivative, hedge)
+                # Negative because selling
+                return -self.criterion.cash(portfolio, target=derivative.payoff())
+
+            mean_price = ensemble_mean(_get_price, n_times=n_times)
 
         return mean_price
diff --git a/pfhedge/nn/modules/loss.py b/pfhedge/nn/modules/loss.py
index fb0aaa53..dc8c5b24 100644
--- a/pfhedge/nn/modules/loss.py
+++ b/pfhedge/nn/modules/loss.py
@@ -4,10 +4,11 @@
 import torch
 from torch import Tensor
 from torch.nn import Module
-from torch.nn import Parameter
+from torch.nn.parameter import Parameter
 
 from pfhedge._utils.bisect import bisect
 from pfhedge._utils.str import _format_float
+from pfhedge._utils.typing import TensorOrScalar
 
 from ..functional import entropic_risk_measure
 from ..functional import exp_utility
@@ -18,24 +19,27 @@
 class HedgeLoss(Module, ABC):
     """Base class for hedging criteria."""
 
-    def forward(self, input: Tensor) -> Tensor:
+    def forward(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
         """Returns the loss of the profit-loss distribution.
 
         This method should be overridden.
 
         Args:
             input (torch.Tensor): The distribution of the profit and loss.
+            target (torch.Tensor or float, default=0): The target portfolio to replicate.
+                Typically, target is the payoff of a derivative.
 
         Shape:
-            - Input: :math:`(N, *)` where
+            - input: :math:`(N, *)` where
               :math:`*` means any number of additional dimensions.
-            - Output: :math:`(*)`
+            - target: :math:`(N, *)`
+            - output: :math:`(*)`
 
         Returns:
             torch.Tensor
         """
 
-    def cash(self, input: Tensor) -> Tensor:
+    def cash(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
         """Returns the cash amount which is as preferable as
         the given profit-loss distribution in terms of the loss.
 
@@ -43,7 +47,7 @@ def cash(self, input: Tensor) -> Tensor:
 
         .. code::
 
-            loss(torch.full_like(pnl, cash)) = loss(pnl)
+            loss(torch.full_like(pl, cash)) = loss(pl)
 
         By default, the output is computed by binary search.
         If analytic form is known, it is recommended to override this method
@@ -51,16 +55,20 @@ def cash(self, input: Tensor) -> Tensor:
 
         Args:
             input (torch.Tensor): The distribution of the profit and loss.
+            target (torch.Tensor or float, default=0): The target portfolio to replicate.
+                Typically, target is the payoff of a derivative.
 
         Shape:
-            - Input: :math:`(N, *)` where
+            - input: :math:`(N, *)` where
               :math:`*` means any number of additional dimensions.
-            - Output: :math:`(*)`
+            - target: :math:`(N, *)`
+            - output: :math:`(*)`
 
         Returns:
             torch.Tensor
         """
-        return bisect(self, self(input), input.min(), input.max())
+        pl = input - target
+        return bisect(self, self(pl), pl.min(), pl.max())
 
 
 class EntropicRiskMeasure(HedgeLoss):
@@ -68,11 +76,11 @@ class EntropicRiskMeasure(HedgeLoss):
     the entropic risk measure.
 
     The entropic risk measure of the profit-loss distribution
-    :math:`\text{pnl}` is given by:
+    :math:`\text{pl}` is given by:
 
     .. math::
-        \text{loss}(\text{pnl}) = \frac{1}{a}
-        \log(- \mathbf{E}[u(\text{pnl})]) \,,
+        \text{loss}(\text{PL}) = \frac{1}{a}
+        \log(- \mathbf{E}[u(\text{PL})]) \,,
         \quad
         u(x) = -\exp(-a x) \,.
 
@@ -86,13 +94,14 @@ class EntropicRiskMeasure(HedgeLoss):
             This parameter should be positive.
 
     Shape:
-        - Input: :math:`(N, *)` where
+        - input: :math:`(N, *)` where
           :math:`*` means any number of additional dimensions.
-        - Output: :math:`(*)`
+        - target: :math:`(N, *)`
+        - output: :math:`(*)`
 
     Examples:
         >>> from pfhedge.nn import EntropicRiskMeasure
-        >>>
+        ...
         >>> loss = EntropicRiskMeasure()
         >>> input = -torch.arange(4.0)
         >>> loss(input)
@@ -101,7 +110,7 @@ class EntropicRiskMeasure(HedgeLoss):
         tensor(-2.0539)
     """
 
-    def __init__(self, a: float = 1.0):
+    def __init__(self, a: float = 1.0) -> None:
         if not a > 0:
             raise ValueError("Risk aversion coefficient should be positive.")
 
@@ -111,20 +120,20 @@ def __init__(self, a: float = 1.0):
     def extra_repr(self) -> str:
         return "a=" + _format_float(self.a) if self.a != 1 else ""
 
-    def forward(self, input: Tensor) -> Tensor:
-        return entropic_risk_measure(input, a=self.a)
+    def forward(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return entropic_risk_measure(input - target, a=self.a)
 
-    def cash(self, input: Tensor) -> Tensor:
-        return -self(input)
+    def cash(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return -self(input - target)
 
 
 class EntropicLoss(HedgeLoss):
     r"""Creates a criterion that measures the expected exponential utility.
 
-    The loss of the profit-loss :math:`\text{pnl}` is given by:
+    The loss of the profit-loss :math:`\text{PL}` is given by:
 
     .. math::
-        \text{loss}(\text{pnl}) = -\mathbf{E}[u(\text{pnl})] \,,
+        \text{loss}(\text{PL}) = -\mathbf{E}[u(\text{PL})] \,,
         \quad
         u(x) = -\exp(-a x) \,.
 
@@ -137,13 +146,14 @@ class EntropicLoss(HedgeLoss):
             the exponential utility.
 
     Shape:
-        - Input: :math:`(N, *)` where
-          :math:`*` means any number of additional dimensions.
-        - Output: :math:`(*)`
+        - input: :math:`(N, *)` where
+            :math:`*` means any number of additional dimensions.
+        - target: :math:`(N, *)`
+        - output: :math:`(*)`
 
     Examples:
         >>> from pfhedge.nn import EntropicLoss
-        >>>
+        ...
         >>> loss = EntropicLoss()
         >>> input = -torch.arange(4.0)
         >>> loss(input)
@@ -152,7 +162,7 @@ class EntropicLoss(HedgeLoss):
         tensor(-2.0539)
     """
 
-    def __init__(self, a: float = 1.0):
+    def __init__(self, a: float = 1.0) -> None:
         if not a > 0:
             raise ValueError("Risk aversion coefficient should be positive.")
 
@@ -162,20 +172,20 @@ def __init__(self, a: float = 1.0):
     def extra_repr(self) -> str:
         return "a=" + _format_float(self.a) if self.a != 1 else ""
 
-    def forward(self, input: Tensor) -> Tensor:
-        return -exp_utility(input, a=self.a).mean(0)
+    def forward(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return -exp_utility(input - target, a=self.a).mean(0)
 
-    def cash(self, input: Tensor) -> Tensor:
-        return -(-exp_utility(input, a=self.a).mean(0)).log() / self.a
+    def cash(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return -(-exp_utility(input - target, a=self.a).mean(0)).log() / self.a
 
 
 class IsoelasticLoss(HedgeLoss):
     r"""Creates a criterion that measures the expected isoelastic utility.
 
-    The loss of the profit-loss :math:`\text{pnl}` is given by:
+    The loss of the profit-loss :math:`\text{PL}` is given by:
 
     .. math::
-        \text{loss}(\text{pnl}) = -\mathbf{E}[u(\text{pnl})] \,,
+        \text{loss}(\text{PL}) = -\mathbf{E}[u(\text{PL})] \,,
         \quad
         u(x) = \begin{cases}
         x^{1 - a} & a \neq 1 \\
@@ -191,13 +201,14 @@ class IsoelasticLoss(HedgeLoss):
             This parameter should satisfy :math:`0 < a \leq 1`.
 
     Shape:
-        - Input: :math:`(N, *)` where
-          :math:`*` means any number of additional dimensions.
-        - Output: :math:`(*)`
+        - input: :math:`(N, *)` where
+            :math:`*` means any number of additional dimensions.
+        - target: :math:`(N, *)`
+        - output: :math:`(*)`
 
     Examples:
         >>> from pfhedge.nn import IsoelasticLoss
-        >>>
+        ...
         >>> loss = IsoelasticLoss(0.5)
         >>> input = torch.arange(1.0, 5.0)
         >>> loss(input)
@@ -206,14 +217,14 @@ class IsoelasticLoss(HedgeLoss):
         tensor(2.3610)
 
         >>> loss = IsoelasticLoss(1.0)
-        >>> pnl = torch.arange(1.0, 5.0)
+        >>> pl = torch.arange(1.0, 5.0)
         >>> loss(input)
         tensor(-0.7945)
         >>> loss.cash(input)
         tensor(2.2134)
     """
 
-    def __init__(self, a: float):
+    def __init__(self, a: float) -> None:
         if not 0 < a <= 1:
             raise ValueError(
                 "Relative risk aversion coefficient should satisfy 0 < a <= 1."
@@ -225,8 +236,8 @@ def __init__(self, a: float):
     def extra_repr(self) -> str:
         return "a=" + _format_float(self.a) if self.a != 1 else ""
 
-    def forward(self, input: Tensor) -> Tensor:
-        return -isoelastic_utility(input, a=self.a).mean(0)
+    def forward(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return -isoelastic_utility(input - target, a=self.a).mean(0)
 
 
 class ExpectedShortfall(HedgeLoss):
@@ -240,13 +251,14 @@ class ExpectedShortfall(HedgeLoss):
             This parameter should satisfy :math:`0 < p \leq 1`.
 
     Shape:
-        - Input: :math:`(N, *)` where
-          :math:`*` means any number of additional dimensions.
-        - Output: :math:`(*)`
+        - input: :math:`(N, *)` where
+            :math:`*` means any number of additional dimensions.
+        - target: :math:`(N, *)`
+        - output: :math:`(*)`
 
     Examples:
         >>> from pfhedge.nn import ExpectedShortfall
-        >>>
+        ...
         >>> loss = ExpectedShortfall(0.5)
         >>> input = -torch.arange(4.0)
         >>> loss(input)
@@ -265,11 +277,11 @@ def __init__(self, p: float = 0.1):
     def extra_repr(self) -> str:
         return str(self.p)
 
-    def forward(self, input: Tensor) -> Tensor:
-        return expected_shortfall(input, p=self.p, dim=0)
+    def forward(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return expected_shortfall(input - target, p=self.p, dim=0)
 
-    def cash(self, input: Tensor) -> Tensor:
-        return -self(input)
+    def cash(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return -self(input - target)
 
 
 class OCE(HedgeLoss):
@@ -278,12 +290,12 @@ class OCE(HedgeLoss):
     The certainty equivalent is given by:
 
     .. math::
-        \\text{loss}(X, w) = w - \\mathrm{E}[u(X + w)]
+        \text{loss}(X, w) = w - \mathrm{E}[u(X + w)]
 
     Minimization of loss gives the optimized certainty equivalent.
 
     .. math::
-        \\rho_u(X) = \\inf_w \\text{loss}(X, w)
+        \rho_u(X) = \inf_w \text{loss}(X, w)
 
     Args:
         utility (callable): Utility function.
@@ -293,13 +305,13 @@ class OCE(HedgeLoss):
 
     Examples:
         >>> from pfhedge.nn.modules.loss import OCE
-        >>>
+        ...
         >>> _ = torch.manual_seed(42)
         >>> m = OCE(lambda x: 1 - (-x).exp())
-        >>> pnl = torch.randn(10)
-        >>> m(pnl)
+        >>> pl = torch.randn(10)
+        >>> m(pl)
         tensor(0.0855, grad_fn=<SubBackward0>)
-        >>> m.cash(pnl)
+        >>> m.cash(pl)
         tensor(-0.0821)
     """
 
@@ -313,5 +325,5 @@ def extra_repr(self) -> str:
         w = float(self.w.item())
         return self.utility.__name__ + ", w=" + _format_float(w)
 
-    def forward(self, input: Tensor) -> Tensor:
-        return self.w - self.utility(input + self.w).mean(0)
+    def forward(self, input: Tensor, target: TensorOrScalar = 0.0) -> Tensor:
+        return self.w - self.utility(input - target + self.w).mean(0)
diff --git a/pfhedge/stochastic/vasicek.py b/pfhedge/stochastic/vasicek.py
index 3b5d63c6..af175bc4 100644
--- a/pfhedge/stochastic/vasicek.py
+++ b/pfhedge/stochastic/vasicek.py
@@ -7,6 +7,8 @@
 
 from pfhedge._utils.typing import TensorOrScalar
 
+from ._utils import cast_state
+
 
 def generate_vasicek(
     n_paths: int,
@@ -72,10 +74,7 @@ def generate_vasicek(
     if init_state is None:
         init_state = (theta,)
 
-    # Accept Union[float, Tensor] as well because making a tuple with a single element
-    # is troublesome
-    if isinstance(init_state, (float, Tensor)):
-        init_state = (torch.as_tensor(init_state),)
+    init_state = cast_state(init_state, dtype, device)
 
     if init_state[0] != 0:
         new_init_state = (init_state[0] - theta,)
@@ -91,10 +90,6 @@ def generate_vasicek(
             device=device,
         )
 
-    # Cast to init_state: Tuple[Tensor, ...] with desired dtype and device
-    init_state = cast(Tuple[Tensor, ...], tuple(map(torch.as_tensor, init_state)))
-    init_state = tuple(map(lambda t: t.to(dtype=dtype, device=device), init_state))
-
     output = torch.empty(*(n_paths, n_steps), dtype=dtype, device=device)
     output[:, 0] = init_state[0]
 
diff --git a/pfhedge/version.py b/pfhedge/version.py
index aa070c2c..5f4bb0b3 100644
--- a/pfhedge/version.py
+++ b/pfhedge/version.py
@@ -1 +1 @@
-__version__ = "0.19.2"
+__version__ = "0.20.0"
diff --git a/pyproject.toml b/pyproject.toml
index e0aa5e82..7416e024 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pfhedge"
-version = "0.19.2"
+version = "0.20.0"
 description = "Deep Hedging in PyTorch"
 authors = ["Shota Imaki <shota.imaki.0801@gmail.com>"]
 license = "MIT"
diff --git a/tests/nn/test_functional.py b/tests/nn/test_functional.py
index f29df3b5..c157a3fb 100644
--- a/tests/nn/test_functional.py
+++ b/tests/nn/test_functional.py
@@ -12,9 +12,9 @@
 from pfhedge.nn.functional import exp_utility
 from pfhedge.nn.functional import expected_shortfall
 from pfhedge.nn.functional import leaky_clamp
+from pfhedge.nn.functional import pl
 from pfhedge.nn.functional import realized_variance
 from pfhedge.nn.functional import realized_volatility
-from pfhedge.nn.functional import terminal_value
 from pfhedge.nn.functional import topp
 from pfhedge.nn.functional import value_at_risk
 
@@ -152,73 +152,75 @@ def test_realized_volatility():
     assert_close(result, expect)
 
 
-def test_terminal_value():
+def test_pl():
     N, T = 10, 20
 
-    # pnl = -payoff if unit = 0
+    # pl = -payoff if unit = 0
     torch.manual_seed(42)
-    spot = torch.randn((N, T)).exp()
-    unit = torch.zeros((N, T))
+    spot = torch.randn((N, 1, T)).exp()
+    unit = torch.zeros((N, 1, T))
     payoff = torch.randn(N)
-    result = terminal_value(spot, unit, payoff=payoff)
+    result = pl(spot, unit, payoff=payoff)
     expect = -payoff
     assert_close(result, expect)
 
     # cost = 0
     torch.manual_seed(42)
-    spot = torch.randn((N, T)).exp()
-    unit = torch.randn((N, T))
-    result = terminal_value(spot, unit)
-    expect = ((spot[..., 1:] - spot[..., :-1]) * unit[..., :-1]).sum(-1)
+    spot = torch.randn((N, 1, T)).exp()
+    unit = torch.randn((N, 1, T))
+    result = pl(spot, unit)
+    expect = ((spot[..., 1:] - spot[..., :-1]) * unit[..., :-1]).sum(-1).squeeze(1)
     assert_close(result, expect)
 
     # diff spot = 0, cost=0 -> value = 0
     torch.manual_seed(42)
-    spot = torch.ones((N, T))
-    unit = torch.randn((N, T))
-    result = terminal_value(spot, unit)
+    spot = torch.ones((N, 1, T))
+    unit = torch.randn((N, 1, T))
+    result = pl(spot, unit)
     expect = torch.zeros(N)
     assert_close(result, expect)
 
     # diff spot = 0, cost > 0 -> value = -cost
     torch.manual_seed(42)
-    spot = torch.ones((N, T))
-    unit = torch.randn((N, T))
-    result = terminal_value(spot, unit, cost=1e-3, deduct_first_cost=False)
-    expect = -1e-3 * ((unit[..., 1:] - unit[..., :-1]).abs() * spot[..., :-1]).sum(-1)
+    spot = torch.ones((N, 1, T))
+    unit = torch.randn((N, 1, T))
+    result = pl(spot, unit, cost=[1e-3], deduct_first_cost=False)
+    expect = -1e-3 * ((unit[..., 1:] - unit[..., :-1]).abs() * spot[..., :-1]).sum(
+        -1
+    ).squeeze(1)
     assert_close(result, expect)
 
     torch.manual_seed(42)
-    spot = torch.ones((N, T))
-    unit = torch.randn((N, T))
-    value0 = terminal_value(spot, unit, cost=1e-3, deduct_first_cost=False)
-    value1 = terminal_value(spot, unit, cost=1e-3, deduct_first_cost=True)
+    spot = torch.ones((N, 1, T))
+    unit = torch.randn((N, 1, T))
+    value0 = pl(spot, unit, cost=[1e-3], deduct_first_cost=False)
+    value1 = pl(spot, unit, cost=[1e-3], deduct_first_cost=True)
     result = value1 - value0
-    expect = -1e-3 * unit[..., 0].abs() * spot[..., 1]
+    expect = -1e-3 * (unit[..., 0].abs() * spot[..., 1]).squeeze(1)
     assert_close(result, expect)
 
 
-def test_terminal_value_unmatched_shape():
-    spot = torch.zeros((10, 20))
-    unit = torch.zeros((10, 20))
+def test_pl_unmatched_shape():
+    spot = torch.zeros((10, 1, 20))
+    unit = torch.zeros((10, 1, 20))
     payoff = torch.zeros(10)
     with pytest.raises(RuntimeError):
-        _ = terminal_value(spot, unit[:-1])
+        _ = pl(spot, unit[:-1])
     with pytest.raises(RuntimeError):
-        _ = terminal_value(spot, unit[:, :-1])
+        _ = pl(spot, unit[:, :-1])
     with pytest.raises(RuntimeError):
-        _ = terminal_value(spot, unit, payoff=payoff[:-1])
+        _ = pl(spot, unit, payoff=payoff[:-1])
 
 
-def test_terminal_value_additional_dim():
+def test_pl_additional_dim():
     N, M, T = 10, 30, 20
 
     # pnl = -payoff if unit = 0
     torch.manual_seed(42)
     spot = torch.randn((N, M, T)).exp()
     unit = torch.zeros((N, M, T))
-    payoff = torch.randn(N, M)
-    result = terminal_value(spot, unit, payoff=payoff)
+    payoff = torch.randn(N)
+    result = pl(spot, unit, payoff=payoff)
     expect = -payoff
     assert_close(result, expect)