Skip to content

Commit 6e2998f

Browse files
George Bisbasgeorgebisbas
George Bisbas
authored and
georgebisbas
committed
examples: Update bs_ivbp
1 parent 8208a53 commit 6e2998f

1 file changed

Lines changed: 60 additions & 58 deletions

File tree

‎examples/finance/bs_ivbp.ipynb‎

Lines changed: 60 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -36,9 +36,11 @@
3636
}
3737
],
3838
"source": [
39-
"from devito import (Eq, Grid, TimeFunction, Operator, solve, Constant,\n",
39+
"from devito import (Eq, Grid, TimeFunction, Operator, solve, Constant, \n",
4040
" SpaceDimension, configuration, centered)\n",
4141
"\n",
42+
"from mpl_toolkits.mplot3d import Axes3D\n",
43+
"from mpl_toolkits.mplot3d.axis3d import Axis\n",
4244
"import matplotlib.pyplot as plt\n",
4345
"import matplotlib as mpl\n",
4446
"\n",
@@ -51,7 +53,7 @@
5153
"\n",
5254
"configuration[\"log-level\"] = 'INFO'\n",
5355
"\n",
54-
"# Constants\n",
56+
"## Constants\n",
5557
"# The strike price of the option\n",
5658
"K = 100.0\n",
5759
"\n",
@@ -67,38 +69,38 @@
6769
"\n",
6870
"# If you want to try some different problems, uncomment these lines\n",
6971
"\n",
70-
"# Example 2\n",
72+
"## Example 2\n",
7173
"# K = 10.0\n",
7274
"# r = 0.1\n",
7375
"# sigma = 0.2\n",
7476
"# smin = 0.0\n",
7577
"# smax = 20.0\n",
7678
"\n",
77-
"# Example 3\n",
79+
"## Example 3\n",
7880
"# K = 100.0\n",
7981
"# r = 0.05\n",
8082
"# sigma = 0.25\n",
8183
"# smin = 50.0\n",
8284
"# smax = 150.0\n",
8385
"\n",
84-
"# Amount of padding to process left/right (hidden from graphs)\n",
86+
"# Amount of padding to proccess left/right (hidden from graphs)\n",
8587
"padding = 10\n",
8688
"smin -= padding\n",
8789
"smax += padding\n",
8890
"\n",
8991
"# Extent calculations\n",
9092
"tmax = 1.0\n",
9193
"dt0 = 0.0005\n",
92-
"ds0 = 1.0\n",
94+
"ds0 = 1.0\n",
9395
"nt = (int)(tmax / dt0) + 1\n",
94-
"ns = int((smax - smin) / ds0) + 1\n",
96+
"ns = int((smax - smin) / ds0) + 1 \n",
9597
"\n",
9698
"shape = (ns, )\n",
97-
"origin = (smin, )\n",
99+
"origin =(smin, )\n",
98100
"spacing = (ds0, )\n",
99101
"extent = int(ds0 * (ns - 1))\n",
100102
"\n",
101-
"print(\"dt,tmax,nt;\", dt0, tmax, nt)\n",
103+
"print(\"dt,tmax,nt;\", dt0,tmax,nt)\n",
102104
"print(\"shape; \", shape)\n",
103105
"print(\"origin; \", origin)\n",
104106
"print(\"spacing; \", spacing)\n",
@@ -162,10 +164,10 @@
162164
"grid = Grid(shape=shape, origin=origin, extent=extent, dimensions=(s, ))\n",
163165
"\n",
164166
"so = 2\n",
165-
"v = TimeFunction(name='v', grid=grid, space_order=so, time_order=1, save=nt)\n",
166-
"v_no_bc = TimeFunction(name='v_no_bc', grid=grid, space_order=so, time_order=1, save=nt)\n",
167+
"v = TimeFunction(name='v', grid=grid, space_order=so, time_order=1, save=nt)\n",
168+
"v_no_bc = TimeFunction(name='v_no_bc', grid=grid, space_order=so, time_order=1, save=nt)\n",
167169
"\n",
168-
"t, s = v.dimensions\n",
170+
"t,s = v.dimensions\n",
169171
"ds = s.spacing\n",
170172
"dt = t.spacing\n",
171173
"\n",
@@ -215,18 +217,18 @@
215217
"outputs": [],
216218
"source": [
217219
"# Equations with Neumann boundary conditions\n",
218-
"eq = [Eq(v[t, extent], v[t, extent-1]+(v[t, extent-1]-v[t, extent-2])),\n",
219-
" Eq(v[t, extent+1], v[t, extent]+(v[t, extent-1]-v[t, extent-2])),\n",
220+
"eq = [Eq(v[t,extent], v[t,extent-1]+(v[t,extent-1]-v[t,extent-2])),\n",
221+
" Eq(v[t,extent+1], v[t,extent]+(v[t,extent-1]-v[t,extent-2])),\n",
220222
" Eq(v.forward, update_centered)]\n",
221223
"eq_no_bc = [Eq(v.forward, update_centered)]\n",
222224
"\n",
223-
"op = Operator(eq, subs=v.grid.spacing_map)\n",
225+
"op = Operator(eq, subs=v.grid.spacing_map)\n",
224226
"op_no_bc = Operator(eq_no_bc, subs=v_no_bc.grid.spacing_map)\n",
225227
"\n",
226228
"# Initial conditions\n",
227229
"\n",
228230
"for i in range(shape[0]):\n",
229-
" v.data[0, i] = max((smin + ds0 * i) - K, 0)\n",
231+
" v.data[0, i] = max((smin + ds0 * i) - K, 0)\n",
230232
" v_no_bc.data[0, i] = max((smin + ds0 * i) - K, 0)"
231233
]
232234
},
@@ -239,8 +241,8 @@
239241
"name": "stderr",
240242
"output_type": "stream",
241243
"text": [
242-
"Operator `Kernel` run in 0.01 s\n",
243-
"Operator `Kernel` run in 0.01 s\n"
244+
"Operator `Kernel` ran in 0.01 s\n",
245+
"Operator `Kernel` ran in 0.01 s\n"
244246
]
245247
},
246248
{
@@ -256,7 +258,7 @@
256258
}
257259
],
258260
"source": [
259-
"# NBVAL_IGNORE_OUTPUT\n",
261+
"#NBVAL_IGNORE_OUTPUT\n",
260262
"\n",
261263
"# Run our operators\n",
262264
"startDevito = timer.time()\n",
@@ -287,27 +289,27 @@
287289
}
288290
],
289291
"source": [
290-
"# NBVAL_IGNORE_OUTPUT\n",
292+
"#NBVAL_IGNORE_OUTPUT\n",
291293
"\n",
292294
"# Get an appropriate ylimit\n",
293295
"slice_smax = v.data[:, int(smax-smin-padding)]\n",
294296
"ymax = max(slice_smax) + 2\n",
295297
"\n",
296298
"# Plot\n",
297299
"s = np.linspace(smin, smax, shape[0])\n",
298-
"plt.figure(figsize=(12, 10), facecolor='w')\n",
300+
"plt.figure(figsize=(12,10), facecolor='w')\n",
299301
"\n",
300302
"time = [1*nt//5, 2*nt//5, 3*nt//5, 4*nt//5, 5*nt//5-1]\n",
301303
"colors = [\"blue\", \"green\", \"gold\", \"darkorange\", \"red\"]\n",
302304
"\n",
303305
"# initial conditions\n",
304-
"plt.plot(s, v_no_bc.data[0, :], '-', color=\"black\", label='initial condition', linewidth=1)\n",
306+
"plt.plot(s, v_no_bc.data[0,:], '-', color=\"black\", label='initial condition', linewidth=1)\n",
305307
"\n",
306308
"for i in range(len(time)):\n",
307-
" plt.plot(s, v_no_bc.data[time[i], :], '-', color=colors[i], label='t='+str(time[i]*dt0), linewidth=1.5)\n",
309+
" plt.plot(s, v_no_bc.data[time[i],:], '-', color=colors[i], label='t='+str(time[i]*dt0), linewidth=1.5)\n",
308310
"\n",
309-
"plt.xlim([smin+padding, smax-padding])\n",
310-
"plt.ylim([0, ymax])\n",
311+
"plt.xlim([smin+padding,smax-padding])\n",
312+
"plt.ylim([0,ymax])\n",
311313
"\n",
312314
"plt.legend(loc=2)\n",
313315
"plt.grid(True)\n",
@@ -341,27 +343,27 @@
341343
}
342344
],
343345
"source": [
344-
"# NBVAL_IGNORE_OUTPUT\n",
346+
"#NBVAL_IGNORE_OUTPUT\n",
345347
"\n",
346348
"# Get an appropriate ylimit\n",
347-
"slice_smax = v.data[:, int(smax-smin-padding)]\n",
349+
"slice_smax = v.data[:,int(smax-smin-padding)]\n",
348350
"ymax = max(slice_smax) + 2\n",
349351
"\n",
350352
"# Plot\n",
351353
"s = np.linspace(smin, smax, shape[0])\n",
352-
"plt.figure(figsize=(12, 10), facecolor='w')\n",
354+
"plt.figure(figsize=(12,10), facecolor='w')\n",
353355
"\n",
354356
"time = [1*nt//5, 2*nt//5, 3*nt//5, 4*nt//5, 5*nt//5-1]\n",
355357
"colors = [\"blue\", \"green\", \"gold\", \"darkorange\", \"red\"]\n",
356358
"\n",
357359
"# initial conditions\n",
358-
"plt.plot(s, v.data[0, :], '-', color=\"black\", label='initial condition', linewidth=1)\n",
360+
"plt.plot(s, v.data[0,:], '-', color=\"black\", label='initial condition', linewidth=1)\n",
359361
"\n",
360362
"for i in range(len(time)):\n",
361-
" plt.plot(s, v.data[time[i], :], '-', color=colors[i], label='t='+str(time[i]*dt0), linewidth=1.5)\n",
363+
" plt.plot(s, v.data[time[i],:], '-', color=colors[i], label='t='+str(time[i]*dt0), linewidth=1.5)\n",
362364
"\n",
363-
"plt.xlim([smin+padding, smax-padding])\n",
364-
"plt.ylim([0, ymax])\n",
365+
"plt.xlim([smin+padding,smax-padding])\n",
366+
"plt.ylim([0,ymax])\n",
365367
"\n",
366368
"plt.legend(loc=2)\n",
367369
"plt.grid(True)\n",
@@ -403,8 +405,9 @@
403405
}
404406
],
405407
"source": [
406-
"# NBVAL_IGNORE_OUTPUT\n",
408+
"#NBVAL_IGNORE_OUTPUT\n",
407409
"\n",
410+
"from mpl_toolkits.mplot3d import Axes3D\n",
408411
"\n",
409412
"# Trim the padding off smin and smax\n",
410413
"trim_data = v.data[:, padding:-padding]\n",
@@ -413,15 +416,15 @@
413416
"tt = np.linspace(0.0, dt0*(nt-1), nt)\n",
414417
"ss = np.linspace(smin+padding, smax-padding, shape[0]-padding*2)\n",
415418
"\n",
416-
"hf = plt.figure(figsize=(12, 12))\n",
419+
"hf = plt.figure(figsize=(12,12))\n",
417420
"ha = plt.axes(projection='3d')\n",
418421
"\n",
419422
"# 45 degree viewpoint\n",
420423
"ha.view_init(elev=25, azim=-45)\n",
421424
"\n",
422-
"ha.set_xlim3d([0.0, 1.0])\n",
423-
"ha.set_ylim3d([smin+padding, smax-padding])\n",
424-
"ha.set_zlim3d([0, ymax])\n",
425+
"ha.set_xlim3d([0.0,1.0])\n",
426+
"ha.set_ylim3d([smin+padding,smax-padding])\n",
427+
"ha.set_zlim3d([0,ymax])\n",
425428
"\n",
426429
"ha.set_xlabel('Time to expiration', labelpad=12, fontsize=16)\n",
427430
"ha.set_ylabel('Stock value', labelpad=12, fontsize=16)\n",
@@ -473,8 +476,8 @@
473476
"name": "stdout",
474477
"output_type": "stream",
475478
"text": [
476-
"devito pde timesteps: 2000, 0.205647s runtime\n",
477-
"call_value_bs timesteps: 5, 3.813932s runtime\n"
479+
"devito pde timesteps: 2000, 0.019737s runtime\n",
480+
"call_value_bs timesteps: 5, 2.596800s runtime\n"
478481
]
479482
},
480483
{
@@ -491,9 +494,9 @@
491494
}
492495
],
493496
"source": [
494-
"# NBVAL_IGNORE_OUTPUT\n",
497+
"#NBVAL_IGNORE_OUTPUT\n",
495498
"\n",
496-
"# Derived formula for Black Scholes call from\n",
499+
"# Derived formula for Black Scholes call from \n",
497500
"# https://aaronschlegel.me/black-scholes-formula-python.html\n",
498501
"def call_value_bs(S, K, T, r, sigma):\n",
499502
" N = Normal('x', 0.0, 1.0)\n",
@@ -504,7 +507,6 @@
504507
" call = (S * cdf(N)(d1) - K * np.exp(-r * T) * cdf(N)(d2))\n",
505508
" return call\n",
506509
"\n",
507-
"\n",
508510
"startBF = timer.time()\n",
509511
"\n",
510512
"# Calculate truth and compare to our solution\n",
@@ -519,21 +521,21 @@
519521
"\n",
520522
"endBF = timer.time()\n",
521523
"\n",
522-
"print(f\"devito pde timesteps: {nt - 1}, {endDevito - startDevito:12.6f}s runtime\")\n",
523-
"print(f\"call_value_bs timesteps: {len(time)}, {endBF - startBF:12.6f}s runtime\")\n",
524-
"\n",
524+
"print(\"devito pde timesteps: %12.6s, %12.6fs runtime\" % (nt-1, endDevito - startDevito))\n",
525+
"print(\"call_value_bs timesteps: %12.6s, %12.6fs runtime\" % (len(time), endBF - startBF))\n",
526+
" \n",
525527
"s2 = np.linspace(smin, smax, shape[0])\n",
526-
"plt.figure(figsize=(12, 10))\n",
528+
"plt.figure(figsize=(12,10))\n",
527529
"\n",
528530
"colors = [\"blue\", \"green\", \"gold\", \"darkorange\", \"red\"]\n",
529-
"plt.plot(s2, v.data[0, :], '-', color=\"black\", label='initial condition', linewidth=1)\n",
531+
"plt.plot(s2, v.data[0,:], '-', color=\"black\", label='initial condition', linewidth=1)\n",
530532
"\n",
531533
"for i in range(len(time)):\n",
532-
" plt.plot(s2, results[i], ':', color=colors[i], label='truth t='+str(time[i]), linewidth=3)\n",
533-
" plt.plot(s2, v.data[int(time[i]*nt), :], '-', color=colors[i], label='pde t='+str(time[i]), linewidth=1)\n",
534+
" plt.plot(s2, results[i], ':', color=colors[i], label='truth t='+str(time[i]), linewidth=3)\n",
535+
" plt.plot(s2, v.data[int(time[i]*nt),:], '-', color=colors[i], label='pde t='+str(time[i]), linewidth=1)\n",
534536
"\n",
535-
"plt.xlim([smin+padding, smax-padding])\n",
536-
"plt.ylim([0, ymax])\n",
537+
"plt.xlim([smin+padding,smax-padding])\n",
538+
"plt.ylim([0,ymax])\n",
537539
"\n",
538540
"plt.legend()\n",
539541
"plt.grid(True)\n",
@@ -560,7 +562,7 @@
560562
{
561563
"data": {
562564
"text/plain": [
563-
"[<matplotlib.lines.Line2D at 0x7f5495584640>]"
565+
"[<matplotlib.lines.Line2D at 0x7fe2669fda50>]"
564566
]
565567
},
566568
"execution_count": 9,
@@ -581,7 +583,7 @@
581583
}
582584
],
583585
"source": [
584-
"# NBVAL_IGNORE_OUTPUT\n",
586+
"#NBVAL_IGNORE_OUTPUT\n",
585587
"\n",
586588
"# Plot the l2 norm of the formula and our solution over time\n",
587589
"t_range = np.linspace(dt0, 1.0, 50)\n",
@@ -592,13 +594,13 @@
592594
" l2 = 0.0\n",
593595
" for x in x_range:\n",
594596
" truth = call_value_bs(x+smin, K, t, r, sigma)\n",
595-
" val = v.data[int(t*(nt-1)), x]\n",
596-
" l2 += (truth - val)**2\n",
597+
" val = v.data[int(t*(nt-1)), x]\n",
598+
" l2 += (truth - val)**2\n",
597599
"\n",
598600
" rms = np.sqrt(np.float64(l2 / len(x_range)))\n",
599601
" vals.append(rms)\n",
600602
"\n",
601-
"plt.figure(figsize=(12, 10))\n",
603+
"plt.figure(figsize=(12,10))\n",
602604
"plt.plot(t_range, np.array(vals))"
603605
]
604606
},
@@ -610,7 +612,7 @@
610612
{
611613
"data": {
612614
"text/plain": [
613-
"0.00581731890853893"
615+
"np.float64(0.005885208362743295)"
614616
]
615617
},
616618
"execution_count": 10,
@@ -619,7 +621,7 @@
619621
}
620622
],
621623
"source": [
622-
"# NBVAL_IGNORE_OUTPUT\n",
624+
"#NBVAL_IGNORE_OUTPUT\n",
623625
"\n",
624626
"np.mean(vals)"
625627
]

0 commit comments

Comments
 (0)