.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/monte_carlo/variance_reduction/plot_05_control_variates.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_monte_carlo_variance_reduction_plot_05_control_variates.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_monte_carlo_variance_reduction_plot_05_control_variates.py:


Lavenberg-Welch control variates
================================

If a quantity :math:`X` simulated alongside the target :math:`Y` has a
known mean, its sampling error tells us about the error in :math:`\bar Y`.
Stephen Lavenberg and Peter Welch (1981) framed the correction as a
regression: estimate :math:`\beta` by least squares and use

.. math:: \hat\mu = \bar Y - \beta^\top(\bar X - \mu_X), \qquad
          \operatorname{Var}\hat\mu \approx (1 - R^2)\frac{s_Y^2}{n}.

The discounted terminal stock price, with mean :math:`S e^{-qT}`, is the
natural control for a call; the closer to the money, the better it
explains the payoff.

.. GENERATED FROM PYTHON SOURCE LINES 19-30

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from financialkit.monte_carlo import control_variate_estimate, terminal_prices, vanilla_payoff
    from financialkit.options import black_scholes_price

    S, T, r, sigma = 100.0, 1.0, 0.05, 0.2
    z = np.random.default_rng(0).standard_normal(10_000)
    S_T = terminal_prices(S, T, r, sigma, z)
    X = np.exp(-r * T) * S_T








.. GENERATED FROM PYTHON SOURCE LINES 31-33

The regression behind the estimator
-----------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 33-48

.. code-block:: Python


    K = 90.0
    Y = np.exp(-r * T) * vanilla_payoff(S_T, K)
    res = control_variate_estimate(Y, X, S)
    print(f"K = {K}: crude {res.crude_estimate:.4f} +/- {res.crude_std_error:.4f}, controlled {res.estimate:.4f} +/- {res.std_error:.4f}")
    print(f"Black-Scholes {black_scholes_price(S, K, T, r, sigma):.4f}, beta {res.beta[0]:.3f}, R^2 {res.r_squared:.3f}")

    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    axes[0].scatter(X[:2000], Y[:2000], s=2, alpha=0.4)
    xs = np.linspace(X.min(), X.max(), 2)
    axes[0].plot(xs, Y.mean() + res.beta[0] * (xs - X.mean()), "k", label=rf"$\beta = {res.beta[0]:.2f}$")
    axes[0].set_xlabel(r"control $e^{-rT} S_T$")
    axes[0].set_ylabel("discounted payoff")
    axes[0].legend()




.. image-sg:: /api/gallery/monte_carlo/variance_reduction/images/sphx_glr_plot_05_control_variates_001.png
   :alt: plot 05 control variates
   :srcset: /api/gallery/monte_carlo/variance_reduction/images/sphx_glr_plot_05_control_variates_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    K = 90.0: crude 16.7734 +/- 0.1731, controlled 16.6782 +/- 0.0418
    Black-Scholes 16.6994, beta 0.835, R^2 0.942

    <matplotlib.legend.Legend object at 0x111f73a10>



.. GENERATED FROM PYTHON SOURCE LINES 49-53

The gain across strikes, with one and two controls
--------------------------------------------------
Adding the discounted at-the-money call as a second control (its
Black-Scholes price is known) helps the out-of-the-money strikes.

.. GENERATED FROM PYTHON SOURCE LINES 53-67

.. code-block:: Python


    strikes = np.linspace(70, 140, 15)
    atm = np.exp(-r * T) * vanilla_payoff(S_T, 100.0)
    controls = np.column_stack([X, atm])
    means = [S, black_scholes_price(S, 100.0, T, r, sigma)]
    one = [1 / (1 - control_variate_estimate(np.exp(-r * T) * vanilla_payoff(S_T, k), X, S).r_squared) for k in strikes]
    two = [1 / (1 - min(control_variate_estimate(np.exp(-r * T) * vanilla_payoff(S_T, k), controls, means).r_squared, 1 - 1e-6)) for k in strikes]
    axes[1].semilogy(strikes, one, "o-", label=r"control $S_T$")
    axes[1].semilogy(strikes, two, "s-", label=r"controls $S_T$ and ATM call")
    axes[1].set_xlabel("strike K")
    axes[1].set_ylabel(r"variance reduction $1/(1 - R^2)$")
    axes[1].legend()
    fig.tight_layout()
    plt.show()








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.080 seconds)


.. _sphx_glr_download_api_gallery_monte_carlo_variance_reduction_plot_05_control_variates.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../../lite/lab/index.html?path=api/gallery/monte_carlo/variance_reduction/plot_05_control_variates.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_05_control_variates.ipynb <plot_05_control_variates.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_05_control_variates.py <plot_05_control_variates.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_05_control_variates.zip <plot_05_control_variates.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_