diff --git a/quantpyml/__init__.py b/quantpyml/__init__.py
index e69de29..4bef480 100644
--- a/quantpyml/__init__.py
+++ b/quantpyml/__init__.py
@@ -0,0 +1,6 @@
+# qol imports
+
+from .returns import brownian_motion, geometric_brownian_motion
+from .equity import EfficientFrontier, BlackScholes
+
+__all__ = ['brownian_motion', 'geometric_brownian_motion', 'EfficientFrontier', 'BlackScholes']
diff --git a/quantpyml/equity/__init__.py b/quantpyml/equity/__init__.py
new file mode 100644
index 0000000..13c8769
--- /dev/null
+++ b/quantpyml/equity/__init__.py
@@ -0,0 +1,6 @@
+# import optimizers
+
+from .efficient_frontier import EfficientFrontier
+from .black_scholes import BlackScholes
+
+__all__ = ['EfficientFrontier', 'BlackScholes']
\ No newline at end of file
diff --git a/quantpyml/optimizers/black_litterman.py b/quantpyml/equity/black_litterman.py
similarity index 100%
rename from quantpyml/optimizers/black_litterman.py
rename to quantpyml/equity/black_litterman.py
diff --git a/quantpyml/equity/black_scholes.py b/quantpyml/equity/black_scholes.py
new file mode 100644
index 0000000..855f8f5
--- /dev/null
+++ b/quantpyml/equity/black_scholes.py
@@ -0,0 +1,49 @@
+
+from scipy.stats import norm
+import numpy as np
+
+class BlackScholes:
+    def __init__(self, price, strike, expiration: "T-days", vol: "annualized constant volatility", rate: "risk free rate" = 0.03, div: "dividend"=0):
+        self.S = price
+        self.K = strike
+        self.T = expiration/365
+        self.r = rate
+        self.sigma = vol
+        self.q = div
+        
+    
+    @staticmethod
+    def N(x):
+        return norm.cdf(x)
+    
+    @property
+    def params(self):
+        return {'S': self.S, 
+                'K': self.K, 
+                'T': self.T, 
+                'r':self.r,
+                'q':self.q,
+                'sigma':self.sigma}
+    
+    def d1(self):
+        return (np.log(self.S/self.K) + (self.r -self.q + self.sigma**2/2)*self.T) / (self.sigma*np.sqrt(self.T))
+    
+    def d2(self):
+        return self.d1() - self.sigma*np.sqrt(self.T)
+    
+    def _call_value(self):
+        return self.S*np.exp(-self.q*self.T)*self.N(self.d1()) - self.K*np.exp(-self.r*self.T) * self.N(self.d2())
+                    
+    def _put_value(self):
+        return self.K*np.exp(-self.r*self.T) * self.N(-self.d2()) - self.S*np.exp(-self.q*self.T)*self.N(-self.d1())
+    
+    def price(self, type_: "call (C), put (P), or both (B)" = 'C'):
+        if type_ == 'C':
+            return self._call_value()
+        if type_ == 'P':
+            return self._put_value() 
+        if type_ == 'B':
+            return  {'call': self._call_value(), 'put': self._put_value()}
+        else:
+            raise ValueError('Unrecognized type')
+
diff --git a/quantpyml/equity/bs.ipynb b/quantpyml/equity/bs.ipynb
new file mode 100644
index 0000000..380535a
--- /dev/null
+++ b/quantpyml/equity/bs.ipynb
@@ -0,0 +1,105 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# reference: https://www.codearmo.com/python-tutorial/options-trading-black-scholes-model\n",
+    "\n",
+    "from scipy.stats import norm\n",
+    "import numpy as np\n",
+    "\n",
+    "class BlackScholes:\n",
+    "    def __init__(self, price, strike, expiration: \"T-days\", vol: \"annualized constant volatility\", rate: \"risk free rate\" = 0.03, div: \"dividend\"=0):\n",
+    "        self.S = price\n",
+    "        self.K = strike\n",
+    "        self.T = expiration/365\n",
+    "        self.r = rate\n",
+    "        self.sigma = vol\n",
+    "        self.q = div\n",
+    "        \n",
+    "    \n",
+    "    @staticmethod\n",
+    "    def N(x):\n",
+    "        return norm.cdf(x)\n",
+    "    \n",
+    "    @property\n",
+    "    def params(self):\n",
+    "        return {'S': self.S, \n",
+    "                'K': self.K, \n",
+    "                'T': self.T, \n",
+    "                'r':self.r,\n",
+    "                'q':self.q,\n",
+    "                'sigma':self.sigma}\n",
+    "    \n",
+    "    def d1(self):\n",
+    "        return (np.log(self.S/self.K) + (self.r -self.q + self.sigma**2/2)*self.T) / (self.sigma*np.sqrt(self.T))\n",
+    "    \n",
+    "    def d2(self):\n",
+    "        return self.d1() - self.sigma*np.sqrt(self.T)\n",
+    "    \n",
+    "    def _call_value(self):\n",
+    "        return self.S*np.exp(-self.q*self.T)*self.N(self.d1()) - self.K*np.exp(-self.r*self.T) * self.N(self.d2())\n",
+    "                    \n",
+    "    def _put_value(self):\n",
+    "        return self.K*np.exp(-self.r*self.T) * self.N(-self.d2()) - self.S*np.exp(-self.q*self.T)*self.N(-self.d1())\n",
+    "    \n",
+    "    def price(self, type_: \"call (C), put (P), or both (B)\" = 'C'):\n",
+    "        if type_ == 'C':\n",
+    "            return self._call_value()\n",
+    "        if type_ == 'P':\n",
+    "            return self._put_value() \n",
+    "        if type_ == 'B':\n",
+    "            return  {'call': self._call_value(), 'put': self._put_value()}\n",
+    "        else:\n",
+    "            raise ValueError('Unrecognized type')\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "{'call': 0.25237702353658875, 'put': 0.17021901090664926}\n"
+     ]
+    }
+   ],
+   "source": [
+    "K = 100\n",
+    "r = 0.1\n",
+    "T = 1\n",
+    "sigma = 0.3\n",
+    "S = 100\n",
+    "print(BlackScholes(S, K, T, r, sigma).price('B'))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "torch",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.16"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/quantpyml/optimizers/efficient_frontier.py b/quantpyml/equity/efficient_frontier.py
similarity index 99%
rename from quantpyml/optimizers/efficient_frontier.py
rename to quantpyml/equity/efficient_frontier.py
index a01a30e..5a231f6 100644
--- a/quantpyml/optimizers/efficient_frontier.py
+++ b/quantpyml/equity/efficient_frontier.py
@@ -1,7 +1,6 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-import sklearn as sk
 import scipy.optimize as spo
 import pandas as pd
 
diff --git a/quantpyml/optimizers/__init__.py b/quantpyml/optimizers/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/quantpyml/returns/__init__.py b/quantpyml/returns/__init__.py
new file mode 100644
index 0000000..0cfb58c
--- /dev/null
+++ b/quantpyml/returns/__init__.py
@@ -0,0 +1,6 @@
+# imports
+
+from .brownian_motion import brownian_motion
+from .geometric_brownian_motion import geometric_brownian_motion
+
+__all__ = ['brownian_motion', 'geometric_brownian_motion']
diff --git a/quantpyml/stochastic/brownian_motion.ipynb b/quantpyml/returns/brownian_motion.ipynb
similarity index 100%
rename from quantpyml/stochastic/brownian_motion.ipynb
rename to quantpyml/returns/brownian_motion.ipynb
diff --git a/quantpyml/returns/brownian_motion.py b/quantpyml/returns/brownian_motion.py
new file mode 100644
index 0000000..75e974c
--- /dev/null
+++ b/quantpyml/returns/brownian_motion.py
@@ -0,0 +1,19 @@
+# implementation of arithmetic brownian motion
+
+import numpy as np
+
+def brownian_motion(nsims=10, ksteps: "number of time steps"=10000, t: "time step" = 0.01, mu: "drift or mean" = 0, sigma: "scaling or standard deviation" = 1):
+    """
+    Weiner process: W(t) ~ N(0,t)
+    Brownian Motion: X(t) = mu*t + sigma*W(t)
+    Standard Brownian Motion: X(t) = W(t)
+    Geometric Brownian Motion: X(t) = exp[(mu-sigma^2/2)*t + sigma*W(t)]
+    """
+    def X():
+        return mu*t + sigma*np.random.normal(loc=0, scale=np.sqrt(t), size=ksteps)
+    x = np.array([i*t for i in range(ksteps)])
+    sims = []
+    for i in range(nsims):
+        y = np.cumsum(X())
+        sims.append(y)
+    return sims, x
diff --git a/quantpyml/returns/geometric_brownian_motion.py b/quantpyml/returns/geometric_brownian_motion.py
new file mode 100644
index 0000000..4ffaf4b
--- /dev/null
+++ b/quantpyml/returns/geometric_brownian_motion.py
@@ -0,0 +1,20 @@
+# implementation of geometric brownian motion
+
+import numpy as np
+
+def geometric_brownian_motion(nsims=10, ksteps: "number of time steps"=1000, t: "time step" = 1, mu: "drift or mean" = 0.0001, sigma: "scaling or standard deviation" = 0.02):
+    """
+    Weiner process: W(t) ~ N(0,t)
+    Brownian Motion: X(t) = X0 + mu*t + sigma*W(t)
+    Standard Brownian Motion: X(t) = X0 + W(t)
+    Geometric Brownian Motion: X(t) = X0 + exp((u-sigma^2/2)*t + sigma*W(t))
+    """
+    def X():
+        return (mu-sigma**2/2)*t + sigma*np.random.normal(loc=0, scale=np.sqrt(t), size=ksteps)
+    x = np.array([i*t for i in range(ksteps)])
+    sims = []
+    for i in range(nsims):
+        y = np.exp(np.cumsum(X()))
+        sims.append(y)
+    return sims, x
+
diff --git a/quantpyml/stochastic/__init__.py b/quantpyml/stochastic/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/test.py b/test.py
new file mode 100644
index 0000000..aa4e0ad
--- /dev/null
+++ b/test.py
@@ -0,0 +1,16 @@
+from quantpyml import brownian_motion
+import matplotlib.pyplot as plt
+import numpy as np
+
+sims, x = brownian_motion()
+# sims, x = geometric_brownian_motion()
+for i in range(len(sims)):
+    plt.plot(x, sims[i])
+# plot a trend line for all the simulations
+plt.plot(x, np.mean(sims, axis=0), color='black', linewidth=3, label='trend')
+plt.xlabel('Time')
+plt.ylabel('Value')
+plt.title('Brownian Motion')
+plt.legend()
+plt.show()
+