The classic textbook theta decay shows that it accelerates until expiry. It is frequently shown with regards to the option value as shown below.
This only holds for ATM options though, because an ITM option will not be worth zero at expiration, and an OTM option will already be worthless some time prior to expiry. Therefore, there will be an inflection point and turnaround somewhere after which theta goes towards zero.
(FD) Finite difference theta (a 1 day bump and reprice theta) is simple to implement for the purpose of figuring out theta and its change over time. You can look here for some details. In Julia, you could use the following code to get a dataframe of theta and its change for various days to expiry, holding all else equal. This code includes the model theta (textbook Black Scholes theta) but the Black Scholes price itself is sufficient if you only use FD theta.
using Plots, Distributions,DataFrames, PlotThemes, PrettyTables, Interact
# define PDF and CDF
N(x) = cdf(Normal(0,1),x)
n(x) = pdf(Normal(0,1),x)
#define Black Scholes with theta (cp = 1 is call, -1 is put)
function BSM(S,K,t,r,d,σ, cp)
d1 = ( log(S/K) + (r - d + 1/2*σ^2)*t ) / (σ*sqrt(t))
d2 = d1 - σ*sqrt(t)
opt = cp*exp(-d*t)S*N(cp*d1) - cp*exp(-r*t)*K*N(cp*d2)
theta = (-(S * exp(-d*t)*n(d1)* σ )/ (2 * sqrt(t)) - r * cp * K * exp(-r*t) * N(cp *d2) +cp*d * S * exp(-d*t)*N(cp *d1))/365
return opt, theta
# define inputs
S,K,r,d,σ,cp_flag = 100,100,0.0,0.02,0.38,1
# compute results
t = 1/365:1/365:42/365
days = 1:1:42
res = BSM.(S,K,reverse(t),r,d,σ,cp_flag)
val = [x for x in res]
theta_bump = append!([val[i+1]-val[i] for i in 1:1:length(val)-1],(cp_flag==1 ? max(S-K,0) : max(K-S,0))+(-val[end]))
theta_change = append!([theta_bump[i+1]-theta_bump[i] for i in 1:1:length(theta_bump)-1], 0)
# output in dataframe
df= DataFrame("Days"=> reverse(days), "TV" => [x for x in res], "Theta" => [x for x in res], "Theta Bumped" => theta_bump, "Theta Acc" => theta_change)
PrettyTables.pretty_table(df, border_crayon = Crayons.crayon"blue", header_crayon = Crayons.crayon"bold green", formatters = ft_printf("%.4f", [2,3,4,5]))
In the dataframe below, TV is the theoretical value of a call, theta is Black Scholes model theta, theta bumped is the finite difference theta, and theta acc is the theta acceleration, which becomes bigger and bigger the closer to expiration in this ATM example (S=K=100).
Since we already have a completely dataframe, we can make an interactive chart using a syntax similar to the one used here. You can clearly see that the shape from the textbook figure above is only valid for approximately ATM options.
Including acceleration, you see the following for an OTM Call:
OTM Put (why theta can be positive here is explained in this answer)
and so forth...
Another observation that may be interesting: Theta crosses the ratio of the price over the number of days left to expiry (P/N) from above at the lowest level of P/N when plotted with time remaining until expiration. This point corresponds graphically to the falling inflection point of theta.
To get this for calendar spreads, you just need to combine two positions.