Brownian motion simulations

Matab code for generating Arithmetic Brownian motion

%This code is written by Tanmay Dutta

function val=ABM_simulation(time,step)

delta=(time/step)^0.5; %This is how we define the delta in BM

p = 0.5; % probability of random walk up or down

val=zeros(1,step);

u=zeros(1,step);

for i=1:numel(u)

u(i)=rand;

%%%%% NOTE ' ' NEAR GREATER AND LESS SIGNS. IF

%% YOU WANT TO REUSE THE CODE REMOVE ' ' IN NEXT LINE

val=delta*(p'<'u)-delta*(p'>'u);

end

for k=2:numel(val)

val(k)=val(k)+val(k-1);

end

val(1)=0.00 % Because B(0) =0 always

end

One simulation result for my code :

(Simulation on T=[0-1] and step = 5000)

(Note that the level goes below zero too, hence it can not be used for stock modeling and we need geometric Brownian motion with drift for that case )

Simulation of different random variables

To simulate a random variable of CDF (note its CDF not PDF) the easiest way is
x = F(inverse) [rand()]
rand() = uniform random variable(matlab and excel has standard functions for this)
Of course this works only when we have some closed form expression for CDF..
Hence will not work for distributions like Gaussian etc.
I will be updating this page for simulating different distributions using matlab

LU factorization for a tridigonal matrix in VBA

Option Base 1
Sub Tridiagonal_solver()
A_Mat = Range("A1:E5")
coeff = Range("H1:H5")
Dim n As Integer
n = UBound(A_Mat, 1)
Dim l_mat() As Double, u_mat() As Double, alpha As Double, beta As Double, y_array() As Double
ReDim y_array(n) As Double
ReDim l_mat(n, n) As Double
ReDim u_mat(n, n) As Double
For i = 1 To n
For j = i To n
l_mat(i, j) = 0
u_mat(i, j) = 0
Next j, i

For i = 1 To n
If (i = 1) Then
u_mat(i, i) = A_Mat(i, i)
End If
l_mat(i, i) = 1

If (i > 1) Then
u_mat(i, i - 1) = A_Mat(i, i - 1)
l_mat(i, i - 1) = A_Mat(i, i - 1) / u_mat(i - 1, i - 1)
u_mat(i, i) = A_Mat(i, i) - l_mat(i, i - 1) * A_Mat(i - 1, i)
End If
Next i
For i = 1 To n
If (i = 1) Then
y_array(i) = coeff(i)
End If

Next i
Range("A8:E12") = l_mat
Range("H8:L12") = u_mat

End Sub

C++ simulation for a puzzle

Sometimes we can find the answers by simulations easier than by logic. (In this case logic is also pretty easy )
The problem : On what day Sat/Sun(For that matters any combination ) will 1st of january fall mostly( Again for that matter any date )
// Purpose : The code try to simulate all the years till 2401 and find which day is on 31st december each year
//Logic:
// Because of the leap year rule we have exacly (400*365+100-3) day and its divisible by 7 completely hence we just
// need to check till 400 years
// Author: Tanmay Dutta
#include
#define start 2002
#define end 2401
using namespace std;
int main()
{
int year,day,count[7];
bool isleap,iscentury;
isleap= false;
iscentury = false;
day=0;
// initializing count to zero
for(int i =0;i<7; i++)
count[i]=0;
enum weekdays{ mon,tue,wed,thr,fri,sat,sun}; // starts with monday = 0
// We know that on 2010 31st dec was on friday
count[day]=1;
for(year = start; year<=end; year++)
{
// check for leap year
if(year%4 ==0)
isleap=true;
// leap year rule does not apply to years with last two digits as 100
if(year%100==0)
isleap = false;
// this rule has an exception at 400 again
if(year%400==0)
isleap=true;
if(isleap)
day=day+2;
else
day=day+1;
if(day==7)
day = 0;
if(day==8)
day=1;
isleap = false;
count[day]+=1;
}
cout<<"Mon"<<'\t'<<"Tue"<<'\t'<<"Wed"<<'\t'<<"Thr"<<'\t'<<"Fri"<<'\t'<<"Sat"<<'\t'<<"Sun"<<'\n';
for(int i =0;i<7; i++)
cout< cout<<'\n';
system("PAUSE");

}

c++ Code for pricing Binary option

The code calculate the bet pricing using the very simplistic Black Scholes model.
The statement for bet is " I want to win X amount if the price of underlying in between Up or Down, While the current price is A"
/*---------------------------------------------------------
#
# Function: Price bets according to BlackScholes framework
# Arguments: Bet amountHigh range, low range, current exchange price,
# time for which bet is valid and implied volatility
# (Can be calculated using current trading options on
# exchange rate )
# Returns: Price for bet in Base currency
#---------------------------------------------------------
*/
#include
#include
#include
using namespace std;
#ifndef Pi
#define Pi 3.141592653589793238462643
#endif
#define interestRate 0.05 // interest rate for the time period
#define maxTries 4 // Number of wrong tries allowed for each entry
double const a0 = 1.26551223, a1 = 1.00002368, a2 = 0.37409196,a3=0.09678418;
double const a4 = 0.18628806, a5 = 0.27886807,a6=1.13520398, a7=1.48851587,a8=0.82215223,a9=0.17087277;
double getCDF(double);
double getArg();

int main()
{
double principle,low_strike,high_strike,time,currentPrice,volatility,bet_price;
cout<<"Enter money to be made out of bet"<<'\n';
principle=getArg();
cout<<"Enter Low range followed by High range and time for bet in days"<<'\n';
low_strike=getArg();high_strike=getArg();time=getArg();
cout<<"Enter the current trading rate and implied volatility(%) for exchange rate(To be determined by fetching data"
"from third party software, currently entered mannually)"<<'\n';
currentPrice=getArg();volatility=getArg();
volatility = volatility/100; // convert to decimal
double years,d_low,d_high;
years= time/365; // covert to years
// Long position
d_low=(log(currentPrice/low_strike) +(interestRate+(volatility*volatility)/2)*years)/(volatility*sqrt(years))-volatility*sqrt(years);
// Short position
d_high=(log(currentPrice/high_strike) +(interestRate+(volatility*volatility)/2)*years)/(volatility*sqrt(years))-volatility*sqrt(years);
bet_price = principle*exp(-1*interestRate*years)*(getCDF(d_low)-getCDF(d_high));
cout< system("PAUSE");
}

double getCDF(double D)
{
double t,erf,z;
z = fabs(D);
z = z/sqrt(2.0);
t = 1/(1+(z*0.5));
erf =t*exp(-z*z-a0+t*(a1+t*(a2+t*(a3+t*(-a4+t*(a5+t*(-a6+t*(a7+
t*(-a8+t*a9)))))))));
erf = erf >= 0.0 ? erf : 2.0-erf;
D = erf/2;
if (D < 0 )
{
D= 1.0 - D;
}
return D;
}
double getArg()
{
double l_arg,exit_c=1;
cin>>l_arg;
while(exit_c {
cin.clear();
cin.ignore(numeric_limits::max(), '\n');
cout << "ERROR: You must enter valid value (no negatives and no alphabets)"<<'\n'<<"<<";
cin >> l_arg;
exit_c+=1;
}
if(exit_c==maxTries)
exit(1);
else
return l_arg;
}

My VBA code for integration (Simpson's Rule)

Sub Simpson()
Sheets("Sheet3").Select
Dim integral As Double, delta As Double, start As Double, last As Double
Dim n As Integer
start = 0
n = 20
last = 1
delta = (last - start) / n
integral = 0
For i = 1 To n / 2
integral = integral + (delta / 3) * (fval(start + (2 * i - 2) * delta) + 4 * fval(start + (2 * i - 1) * delta) + fval(start + 2 * i * delta))


Next
Cells(1, 1) = integral
End Sub

My VBA code for integration (Trapezoidal as approximation of area)

Sub Trapeziodal()
Sheets("Sheet2").Select
Dim integral As Double, delta As Double, start As Double, last As Double
Dim n As Integer
start = 0
n = 10
last = 1
delta = (last - start) / n
integral = 0
For i = 1 To n
integral = integral + 0.5 * delta * (fval(start + (i - 1) * delta) + fval(start + (i * delta)))


Next
Cells(1, 1) = integral
End Sub

My VBA code for integration (Rectangle as approximation of area)

Sub rectangle()
Sheets("Sheet1").Select
Dim integral As Double, delta As Double, start As Double, last As Double
Dim n As Integer
start = 0
n = 10
last = 1
delta = (last - start) / n
integral = 0
For i = 1 To n
integral = integral + delta * fval(start + (i - 1) * delta)

Next
Cells(1, 1) = integral
End Sub