.. 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_06_moment_matching.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_06_moment_matching.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_06_moment_matching.py:


Barraquand moment matching
==========================

A finite sample of normals never has mean exactly 0 and variance exactly
1, and the payoff inherits those sampling errors. Jérôme Barraquand
(1995), pricing multi-asset options, simply rescaled the sample,

.. math:: \tilde z_i = \frac{z_i - \bar z}{s_z},

so the first two moments are exact ("quadratic resampling"). The draws
become dependent, so error bars come from independent replications.

.. GENERATED FROM PYTHON SOURCE LINES 16-31

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from financialkit.monte_carlo import moment_match, moment_matching_estimate, sample_mean, terminal_prices, vanilla_payoff
    from financialkit.monte_carlo.visualizers import plot_estimator_comparison
    from financialkit.options import black_scholes_price

    S, K, T, r, sigma = 100.0, 100.0, 1.0, 0.05, 0.2
    exact = black_scholes_price(S, K, T, r, sigma)


    def call(z):
        return np.exp(-r * T) * vanilla_payoff(terminal_prices(S, T, r, sigma, z), K)









.. GENERATED FROM PYTHON SOURCE LINES 32-34

Matched samples have exact moments
----------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 34-40

.. code-block:: Python


    z = np.random.default_rng(0).standard_normal(500)
    w = moment_match(z)
    print(f"raw:     mean {z.mean():+.4f}, std {z.std():.4f}")
    print(f"matched: mean {w.mean():+.4f}, std {w.std():.4f}")





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

 .. code-block:: none

    raw:     mean -0.0269, std 1.0136
    matched: mean -0.0000, std 1.0000




.. GENERATED FROM PYTHON SOURCE LINES 41-43

Twenty replications of 2000 draws each
--------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 43-53

.. code-block:: Python


    n, reps = 2000, 20
    gen = np.random.default_rng(1)
    crude = sample_mean([call(gen.standard_normal(n)).mean() for _ in range(reps)], method="crude", n_evaluations=n * reps)
    matched = moment_matching_estimate(call, n, n_replications=reps, rng=2)
    print(f"crude {crude.estimate:.4f} +/- {crude.std_error:.4f}; matched {matched.estimate:.4f} +/- {matched.std_error:.4f}; exact {exact:.4f}")

    fig, axes = plt.subplots(1, 2, figsize=(11, 4))
    plot_estimator_comparison([crude, matched], exact=exact, ax=axes[0])




.. image-sg:: /api/gallery/monte_carlo/variance_reduction/images/sphx_glr_plot_06_moment_matching_001.png
   :alt: plot 06 moment matching
   :srcset: /api/gallery/monte_carlo/variance_reduction/images/sphx_glr_plot_06_moment_matching_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    crude 10.2772 +/- 0.0496; matched 10.4502 +/- 0.0072; exact 10.4506

    <Axes: ylabel='estimate (95% CI); variance reduction'>



.. GENERATED FROM PYTHON SOURCE LINES 54-56

The gain is largest for small samples
-------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 56-68

.. code-block:: Python


    ns = np.array([100, 300, 1000, 3000, 10000])
    ratio = []
    for k in ns:
        g = np.random.default_rng(int(k))
        c = sample_mean([call(g.standard_normal(k)).mean() for _ in range(reps)], n_evaluations=k * reps)
        ratio.append(moment_matching_estimate(call, int(k), n_replications=reps, rng=int(k) + 1).variance_ratio(c))
    axes[1].semilogx(ns, ratio, "o-")
    axes[1].set_xlabel("draws per replication")
    axes[1].set_ylabel("variance reduction factor")
    fig.tight_layout()
    plt.show()








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

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


.. _sphx_glr_download_api_gallery_monte_carlo_variance_reduction_plot_06_moment_matching.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_06_moment_matching.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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