#!/usr/bin/perl
require 5.008;
our $VERSION = "0.01"; # Time-stamp: <2017-02-05T12:55:29Z>

use strict;
use warnings;
use utf8;
use Math::Trig;
use Math::Gauss;

our $TRIALS = 1000;		# 試行数
our $VOLATILITY = 0.1;		# ボラティリティ(σ または θ)
our $INTEREST = 0.06;		# 無危険利子率
our $INITIAL_PRICE = 18000;	# 初期価格
our $STRIKE_PRICE = 19000;	# 行使価格
our $TERM = 5;			# 期間
our $STEP = 10;			# 1期間の分割数
our $MODE = 0;

use Getopt::Long;
use Pod::Usage;
$ENV{"PERLDOC"} = "" if ! exists $ENV{"PERLDOC"};
$ENV{"PERLDOC"} .= " " if $ENV{"PERLDOC"} ne "";
$ENV{"PERLDOC"} .= "-wcenter:'Probability Model'";

Getopt::Long::Configure("bundling", "auto_version");
GetOptions(
	   "trial=i" => \$TRIALS,
	   "volatility=f" => \$VOLATILITY,
	   "interest=f" => \$INTEREST,
	   "initial=f" => \$INITIAL_PRICE,
	   "strike=f" => \$STRIKE_PRICE,
	   "term=i" => \$TERM,
	   "step=i" => \$STEP,
	   "normal" => sub {$MODE = 1;},
	   "cauchy" => sub {$MODE = 2;},

	   "man" => sub {pod2usage(-verbose => 2)},
	   "h|?" => sub {pod2usage(-verbose => 0, -output=>\*STDOUT, 
				   -exitval => 1)},
	   "help" => sub {pod2usage(1)},
	   ) or pod2usage(-verbose => 0);

my $U = ($INTEREST / $STEP) + ($VOLATILITY / sqrt($STEP));
my $D = ($INTEREST / $STEP) - ($VOLATILITY / sqrt($STEP));

MAIN_LOOP:
{
  our $accume = 0;

  for (my $i = 0; $i < $TRIALS; $i++) {
    my $price = $INITIAL_PRICE;
    for (my $i = 0; $i < $TERM * $STEP; $i++) {
      my $x;
      if ($MODE == 1) {
	my $r1 = sqrt(-2 * log(rand(1))) * cos(2 * pi * rand(1));
	$x = 1 + ($INTEREST / $STEP) + ($VOLATILITY / sqrt($STEP)) * $r1;
      } elsif ($MODE == 2) {
	my $r1 = tan(pi * (rand(1) - 0.5));
	$x = 1 + ($INTEREST / $STEP) + ($VOLATILITY / sqrt($STEP)) * $r1;
      } else {
	$x = (rand(1) < 0.5)? (1 + $U) : (1 + $D);
      }
      $price *= $x;
    }
    my $profit = ($price > $STRIKE_PRICE)? ($price - $STRIKE_PRICE) : 0;
    my $discounted = $profit / ((1 + ($INTEREST / $STEP)) ** ($TERM * $STEP));
#    print "d : $discounted $profit $price\n";
    $accume += $discounted;
  }
  printf("result = %g\n", $accume / $TRIALS);

  my $tmp1 = log($INITIAL_PRICE / $STRIKE_PRICE) + $TERM * $INTEREST;
  my $tmp2 = $TERM * (($VOLATILITY ** 2) / 2);
  my $tmp3 = $VOLATILITY * sqrt($TERM);
  my $tmp4 = Math::Gauss::cdf(($tmp1 + $tmp2) / $tmp3);
  my $tmp5 = Math::Gauss::cdf(($tmp1 - $tmp2) / $tmp3);
  my $tmp = $INITIAL_PRICE * $tmp4 
    - $STRIKE_PRICE * exp(- $INTEREST * $TERM) * $tmp5;
  printf("estimated = %g\n" , $tmp);
}

__END__

=pod

=encoding utf8

=head1	NAME

test_bs_1.pl - 2項モデルによるブラック＝ショールズ式の近似

=head1	SYNOPSIS

perl B<test_bs_1.pl> --trial=TRIALS --volatility=VOLATILITY [--normal|--cauchy]

=head1	Options

=over 8

=item B<--help>

show help message about options.

=item B<--man>

show man page.

=item B<--version>

show version infomation.

=item B<--trial>

試行数。

=item B<--normal>

2項モデルで単純な確率1/2の上下の変動の替わりに正規乱数を使ってみる。

=item B<--cauchy>

2項モデルで単純な確率1/2の上下の変動の替わりにコーシー分布の乱数を使ってみる。

=item B<--volatility>

ボラティリティ σ (または θ) を指定する。

=back

=head1	DESCRIPTION

まず、2項モデルをモンテカルロ法で実行し、ブラック＝ショールズ式で近似できることを確かめる。

次に、2項モデルの上下二択の変わりに、正規分布やコーシー分布の乱数を使ってみる。すると、正規分布だと元のモデルとオプション価格が変わらないのに対し、コーシー分布だと大きく変わることがわかる。デフォルトの数値例で、σ = 0.1 に対応するようなコーシー分布のθ(このプログラムでは σ で指定)は、0.00005 ぐらいでないとダメなようだ。

ただし、--step を少なくすれば、コーシー分布も落ち着いてくる。そこが正規分布を使ったときとの違い。つまり、この実験で、コーシー分布をブラック＝ショールズ・モデル風に使うのはうまくいっていないということだろう。

=head1	AUTHOR

JRF

=head1	LICENSE

Public Domain.

=cut