<a href="https://colab.research.google.com/github/nickwotton/MQP2019/blob/master/Nick/BSM_NN_v02.ipynb" target="_parent"></a>

# Attempt to Replicate the Black Scholes Model Using a Neural Network



```python
import torch
import torch.nn as nn
import numpy as np
import torch.nn.functional as F
import matplotlib.pyplot as plt
import scipy.stats as ss
```

## Define the Function
Here we define our function, the Black Scholes Model (BSM). First, we must initialize the option class, then the Geometric Brownian Motion Class, and finally the BSM class.
Then we test the equation with a test value of 2.


```python
'''=========
option class init
=========='''
class VanillaOption:
    def __init__(
        self,
        otype = 1, # 1: 'call'
                  # -1: 'put'
        strike = 110.,
        maturity = 1.,
        market_price = 10.):
      self.otype = otype               # Put or Call
      self.strike = strike             # Strike K
      self.maturity = maturity         # Maturity T
      self.market_price = market_price #this will be used for calibration
      
        
    def payoff(self, s): #s: excercise price
      otype = self.otype
      k = self.strike
      maturity = self.maturity
      return np.max([0, (s - k)*otype])

```


```python
'''============
Gbm class
============='''

class Gbm:
    def __init__(self,
                 init_state = 100.,
                 drift_ratio = .0475,
                 vol_ratio = .2
                ):
        self.init_state = init_state
        self.drift_ratio = drift_ratio
        self.vol_ratio = vol_ratio
```


```python
'''========
Black-Scholes-Merton formula. 
=========='''

def bsm_price(self, vanilla_option):
    s0 = self.init_state
    sigma = self.vol_ratio
    r = self.drift_ratio
    
    otype = vanilla_option.otype
    k = vanilla_option.strike
    maturity = vanilla_option.maturity
    
    d1 = 1/(sigma*np.sqrt(maturity))*(np.log(s0/k) + (r + np.power(sigma,2)/2)*(maturity)) 
    d2 = 1/(sigma*np.sqrt(maturity))*(np.log(s0/k) + (r - np.power(sigma,2)/2)*(maturity)) 
    return (otype * s0 * ss.norm.cdf(otype * d1) #line break needs parenthesis
            - otype * np.exp(-r * maturity) * k * ss.norm.cdf(otype * d2))

Gbm.bsm_price = bsm_price
```


```python
'''=======
Get BSM prices given an option and a Tensor
======='''

def prices_bsm(self, vanilla_option, data):
  
  # Create the list
  a = []

  # Get tensor size
  sizeT = list(data.size())[0]

  for i in range(sizeT):
    self.init_state = data[i].item()
    callPrice = gbm1.bsm_price(vanilla_option)
    a.append(callPrice)
  
  # Convert to array then tensor
  arrayOut = np.array(a)
  outputData = torch.from_numpy(arrayOut)
  return outputData

Gbm.prices_bsm = prices_bsm
```


```python
def f(s):
  gbm = Gbm(init_state=s)
  option = VanillaOption(strike=10)
  return gbm.bsm_price(option)

batch_size = 21
x_list = np.linspace(8, 11, batch_size)
y_list = np.array([f(x) for x in x_list])
#plt.plot(x_list, y_list)
```

## Create Model
Next, we create the neural network model. This is done first by setting the inner and outer dimensions with variables. Next we code the model and vary the internal dimensions to attempt to improve the model. At this level, this is essentially a simple linear algebra exercise:

If we have input $x$, internal parameters $a,b$, and solution $f(x)$ then in the one-dimensional case we have:

\begin{equation}
  \left(
    a_{1}x+b_{2}
  \right)
  a_{2} + b_{2}
  = f(x)
\end{equation}
  
However, we want to get a better estimate for the true equation. So we increase the interior dimension which corresponds to the number of neurons inside the network. For example, we raised the inner dimension to 3. In matrix form we have:

\begin{equation}
\left(
  \begin{bmatrix} x \end{bmatrix} 
  \begin{bmatrix} a_{1} & a_{2} & a_{3} \end{bmatrix} 
  + 
  \begin{bmatrix} b_{1} & b_{2} & b_{3} \end{bmatrix}
\right)
 \begin{bmatrix} a_{4} \\ a_{5} \\ a_{6} \end{bmatrix}
 +
 \begin{bmatrix} b_{4} \\ \end{bmatrix}
 =
 \begin{bmatrix} f(x) \end{bmatrix}
\end{equation}

Graphically, we can render this second neural network as:


**ADD DESCRIPTION**


```python
#model
#nn.Linear
H1 = 50; H2 = 11 #number of hidden layer
model = nn.Sequential(
    nn.Linear(1, H1), 
    nn.Sigmoid(),
    nn.Linear(H1, H2),
    nn.Sigmoid(),
    nn.Linear(H2,2),
    nn.Sigmoid(),
    nn.Linear(2,1)
)
```

Here we define the Loss function as the Mean Squared Error(MSE). 

Note that by doing so, we are essentially 'cheating' the system. In most applications, we would not know the function $f$ so we would be unable to find the MSE.


```python
#loss function 
criterion = nn.MSELoss()
```

Next we choose a learning rate and a method for learning. The learning rate is the percent of the data that is accepted in each iteration. The Methods we tried were SGD and Adam. 


```python
#optimizer
learning_rate = 0.1
optimizer = torch.optim.SGD(model.parameters(), lr=learning_rate, momentum=0.7) 
```


```python
batch_size = np.size(x_list)
x_train0 = torch.from_numpy(x_list).reshape(batch_size,1).float()
y_train0 = torch.from_numpy(y_list).reshape(batch_size,1).float()
```

## Train the Model

First we create the training data. This is a batch of random points that we pass through the BSM.


```python
# Create Training Data

#batch_size = 1000

# Create a list
#a = []

# Generate randomized data
#x_train = torch.randn(batch_size, 1)
# Modify the random data to be a stock price
#x_train0 = np.abs(x_train*100)

# Create static training data for conistency
#x_train = torch.tensor([[100], [105], [110], [115], [120]])

# Initialize GBM
#gbm1 = Gbm()
#option = VanillaOption()

# Call the prices_bsm to get y_train
#y_train0 = gbm1.prices_bsm(option, x_train)

#print(x_train[1:10])
#print(y_train[1:10])

```

Once we have the training data, we pass this collection of inputs and solutions into the model. With each iteration we calculate the loss and attempt to optimize the model to further reduce the loss.



```python
# Normalization

def linear_transform(xx, l = 0, u= 1):
  M = torch.max(xx)
  m = torch.min(xx)
  return (u-l)/(M-m)*(xx-m)+l, m, M, l, u
x_train, x_m, x_M, x_l, x_u = linear_transform(x_train0, -1, 1)
y_train, y_m, y_M, y_l, y_u = linear_transform(y_train0, 0, 1)
```


```python
# Train the model

num_epochs = 1000

for epoch in range(num_epochs):

    # Forward pass
    outputs = model(x_train.float())
    loss = criterion(outputs, y_train.float())
    
    # Backward and optimize
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
    
    if epoch == 0 or (epoch+1) % 50 == 0:
        print ('Epoch [{}/{}], Loss: {:.4f}'.format(epoch+1, 
                                                    num_epochs, loss.item()))
        print(outputs[1:10])
       
      

```

    Epoch [1/1000], Loss: 0.1023
    tensor([[0.3266],
            [0.3266],
            [0.3265],
            [0.3264],
            [0.3264],
            [0.3263],
            [0.3263],
            [0.3262],
            [0.3261]], grad_fn=<SliceBackward>)
    Epoch [50/1000], Loss: 0.0950
    tensor([[0.4094],
            [0.4094],
            [0.4094],
            [0.4094],
            [0.4094],
            [0.4094],
            [0.4095],
            [0.4095],
            [0.4095]], grad_fn=<SliceBackward>)
    Epoch [100/1000], Loss: 0.0946
    tensor([[0.4082],
            [0.4084],
            [0.4085],
            [0.4086],
            [0.4087],
            [0.4088],
            [0.4090],
            [0.4091],
            [0.4092]], grad_fn=<SliceBackward>)
    Epoch [150/1000], Loss: 0.0941
    tensor([[0.4068],
            [0.4071],
            [0.4073],
            [0.4076],
            [0.4078],
            [0.4081],
            [0.4084],
            [0.4086],
            [0.4089]], grad_fn=<SliceBackward>)
    Epoch [200/1000], Loss: 0.0933
    tensor([[0.4047],
            [0.4052],
            [0.4056],
            [0.4061],
            [0.4065],
            [0.4070],
            [0.4075],
            [0.4079],
            [0.4084]], grad_fn=<SliceBackward>)
    Epoch [250/1000], Loss: 0.0920
    tensor([[0.4011],
            [0.4019],
            [0.4027],
            [0.4035],
            [0.4043],
            [0.4051],
            [0.4060],
            [0.4068],
            [0.4077]], grad_fn=<SliceBackward>)
    Epoch [300/1000], Loss: 0.0894
    tensor([[0.3942],
            [0.3956],
            [0.3971],
            [0.3986],
            [0.4002],
            [0.4017],
            [0.4033],
            [0.4049],
            [0.4065]], grad_fn=<SliceBackward>)
    Epoch [350/1000], Loss: 0.0836
    tensor([[0.3789],
            [0.3819],
            [0.3849],
            [0.3880],
            [0.3911],
            [0.3943],
            [0.3975],
            [0.4008],
            [0.4040]], grad_fn=<SliceBackward>)
    Epoch [400/1000], Loss: 0.0697
    tensor([[0.3400],
            [0.3468],
            [0.3538],
            [0.3609],
            [0.3682],
            [0.3757],
            [0.3832],
            [0.3908],
            [0.3985]], grad_fn=<SliceBackward>)
    Epoch [450/1000], Loss: 0.0384
    tensor([[0.2362],
            [0.2529],
            [0.2703],
            [0.2883],
            [0.3069],
            [0.3260],
            [0.3454],
            [0.3652],
            [0.3852]], grad_fn=<SliceBackward>)
    Epoch [500/1000], Loss: 0.0077
    tensor([[0.0706],
            [0.1010],
            [0.1335],
            [0.1679],
            [0.2042],
            [0.2422],
            [0.2816],
            [0.3222],
            [0.3637]], grad_fn=<SliceBackward>)
    Epoch [550/1000], Loss: 0.0020
    tensor([[-0.0010],
            [ 0.0325],
            [ 0.0689],
            [ 0.1083],
            [ 0.1507],
            [ 0.1958],
            [ 0.2436],
            [ 0.2936],
            [ 0.3454]], grad_fn=<SliceBackward>)
    Epoch [600/1000], Loss: 0.0011
    tensor([[-0.0026],
            [ 0.0284],
            [ 0.0625],
            [ 0.0998],
            [ 0.1403],
            [ 0.1841],
            [ 0.2310],
            [ 0.2808],
            [ 0.3332]], grad_fn=<SliceBackward>)
    Epoch [650/1000], Loss: 0.0007
    tensor([[0.0047],
            [0.0333],
            [0.0650],
            [0.0999],
            [0.1382],
            [0.1799],
            [0.2251],
            [0.2735],
            [0.3250]], grad_fn=<SliceBackward>)
    Epoch [700/1000], Loss: 0.0005
    tensor([[0.0103],
            [0.0373],
            [0.0673],
            [0.1005],
            [0.1372],
            [0.1773],
            [0.2211],
            [0.2684],
            [0.3191]], grad_fn=<SliceBackward>)
    Epoch [750/1000], Loss: 0.0003
    tensor([[0.0144],
            [0.0402],
            [0.0690],
            [0.1010],
            [0.1364],
            [0.1755],
            [0.2182],
            [0.2647],
            [0.3148]], grad_fn=<SliceBackward>)
    Epoch [800/1000], Loss: 0.0003
    tensor([[0.0176],
            [0.0425],
            [0.0704],
            [0.1014],
            [0.1359],
            [0.1741],
            [0.2161],
            [0.2619],
            [0.3114]], grad_fn=<SliceBackward>)
    Epoch [850/1000], Loss: 0.0002
    tensor([[0.0201],
            [0.0443],
            [0.0715],
            [0.1018],
            [0.1356],
            [0.1731],
            [0.2144],
            [0.2597],
            [0.3088]], grad_fn=<SliceBackward>)
    Epoch [900/1000], Loss: 0.0002
    tensor([[0.0220],
            [0.0458],
            [0.0724],
            [0.1022],
            [0.1354],
            [0.1724],
            [0.2132],
            [0.2580],
            [0.3068]], grad_fn=<SliceBackward>)
    Epoch [950/1000], Loss: 0.0001
    tensor([[0.0237],
            [0.0470],
            [0.0731],
            [0.1025],
            [0.1353],
            [0.1718],
            [0.2122],
            [0.2567],
            [0.3052]], grad_fn=<SliceBackward>)
    Epoch [1000/1000], Loss: 0.0001
    tensor([[0.0249],
            [0.0479],
            [0.0737],
            [0.1027],
            [0.1352],
            [0.1713],
            [0.2114],
            [0.2556],
            [0.3039]], grad_fn=<SliceBackward>)


## Testing the Model



```python
def learnedfun(x):
  out = (1-(-1))/(x_M-x_m)*(x-x_m)+(-1.)
  out = model(out)
  out = (y_M- y_m)*out+y_m
  return out

y_pred = learnedfun(x_train0)
```


```python
# Test with training data
plt.scatter(x_train0.detach().numpy(), y_train0.detach().numpy(), label='true')
plt.plot(x_train0.detach().numpy(), y_pred.detach().numpy(), label='pred')

plt.legend()
plt.show()
```


```python
# Test the Model with Random data

# Generate random data
x_0 = torch.randn(50,1)

# Get BSM Prices as determined by the formula
y_0 = gbm1.prices_bsm(option, x_)

# Normalize
x_, x_m, x_M, x_l, x_u = linear_transform(x_0)
y_, y_m, y_M, y_l, y_u = linear_transform(x_0)

# Plot x_ versus formula prices
plt.scatter(x_0.detach().numpy(), y_0.detach().numpy(), label='true')

# Get BSM Prices as determined by the model
y_pred = learnedfun(x_0)

# Plot x_ versus the model prices
plt.scatter(x_0.detach().numpy(), y_pred.detach().numpy(), label='pred')

plt.legend()

```


```python

```

a954d893079cb81c62c8cb9950b462c337f5cae7
ipynb
Jupyter Notebook
Nick/BSM_NN_v02.ipynb
nickwotton/MQP2019
8d7e7db5d47d605a605aeb394cf5b8af2b13e7fb
MIT
2019-09-13T16:27:14.000Z
2019-09-13T16:27:14.000Z
Nick/BSM_NN_v02.ipynb
nickwotton/MQP2019
8d7e7db5d47d605a605aeb394cf5b8af2b13e7fb
MIT
2019-09-02T17:49:28.000Z
2020-04-06T13:03:43.000Z
Nick/BSM_NN_v02.ipynb
nickwotton/MQP2019
8d7e7db5d47d605a605aeb394cf5b8af2b13e7fb
MIT
2019-09-11T01:51:56.000Z
2019-12-03T20:13:16.000Z
Qwen/Qwen-72B
1. YES
2. YES