The OP is no longer active, but a lot of people looked at this question, which is based on a misunderstanding of what theta is. If you type the formula in Wikipedia into excel, you should see that it matches your data provider. Technically, SPY options are American, but given the option at hand, there is no difference to a European option anyways. In general, the difference, especially in theta will be very small.
Most users here will not see closed questions on quant.stackexchange, but if you do, you can look here. The screenshot with the values looks like this:
Instead of Excel, I will use Julia because charting will be easier and interactive. In any case, you can almost copy paste the code for the Greeks and option value from the wikipedia link provided in the question. As mentioned there, theta is expressed in value per year and usually divided by the number of days in a year.
using Interact,Plots, Distributions,DataTheta Frames, PlotThemes, Dates, PrettyTables, Interact
N(x) = cdf(Normal(0,1),x)
n(x) = pdf(Normal(0,1),x)
Calculate Black-Scholes European call option price
d1 = ( log(S/K) + (r - d + 1/2*σ^2)*t ) / (σ*sqrt(t))
d2 = d1 - σ*sqrt(t)
p = -exp(-d*t)S*N(-d1) + exp(-r*t)*K*N(-d2)
delta_p = -exp(-d*t)*N(-d1)
gamma_p = exp(-d*t)*n(d1) / (S*σ *sqrt(t))
theta_p = (-(S * exp(-d*t)*n(d1)* σ )/ (2 * sqrt(t)) + r * K * exp(-r*t) * N(-d2) - d * S * exp(-d*t)*N(-d1))/365
rho_p = -( K*t * exp(-d*t) * N(-d2))*0.01
vega_p = S * exp(-d*t)*n(d1) * sqrt(t)*0.01
return delta_p, gamma_p, theta_p, vega_p, rho_p, p
df= DataFrame("Days"=> reverse(days), "Delta" => [x for x in res], "Gamma" => [x for x in res], "Theta" => [x for x in res], "Vega" => [x for x in res], "Rho" => [x for x in res], "Theoretical value" => [x for x in res], "Theta Bumped" => theta_bump, "P/N" => p_n)
hl_1 = Highlighter((data,i,j) -> data[i,1] == 42, crayon"bg:dark_gray white bold")
PrettyTables.pretty_table(df, border_crayon = Crayons.crayon"blue", header_crayon = Crayons.crayon"bold green", formatters = ft_printf("%.4f", [2,3,4,5,6]), highlighters = (hl_1))
The dataframe below is the result, with the first line being the option priced in the screenshot above. As you can see, it unsurprisingly (it is just Black Scholes) matches exactly.
The other lines shorten time to expiry by one day at a time, keeping all else equal. Theta bumped is the finite difference computation of theta (shorten by one day, look at the difference as shown here. You simply cannot multiply theta by the number of days because theta changes with time (t is part of the formula).
With regards to P/N, that also does not work because theta shows the change in price for one day, not some average per day, spread equally over the options life.
The problem with the chart, and the interpretation that
one would expect that at any given point in time, the theta of an
option would generally be less than its current price, divided by the
days remaining; the reason being, it's going to continue to decay
is that looking blindly at a single figure tells you little about option pricing. If you think a bit about the figure, you will quickly realize that it can only be for an ATM option 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 somewhere after which theta goes towards zero (as in the dataframe above). Since we already have a completely dataframe, we can make an interactive chart using a syntax similar to the one used here.
So is there any relation to P/N at all? Yes! Theta crosses 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. Where this occurs depends a lot on moneyness and IV, similar to any other value in option pricing. The vertical blue line is the 42 days of the option in the example of the OP.
The Natenberg example mentioned in another answer is from futures, which does not apply to ETFs (SPY). Natenberg however has a similar argument for stock options on P. 109 which goes like this
Is it ever possible for an option to have a positive theta such that if nothing
changes, the option will be worth more tomorrow than it is today? In fact, this
can happen because of the depressing effect of interest rates. Consider a 60 call
on an underlying contract that is currently trading at 100. How much might
this call be worth if we know that at expiration the underlying contract will still
be at 100? At expiration, the option will be worth 40, its intrinsic value. However, if the option is subject to stock-type settlement, today it will only be worth
the present value of 40, perhaps 39. If the underlying price remains at 100, as
time passes, the value of the option must rise from 39 (its value today) to 40 (its
intrinsic value at expiration). The option in effect has negative time value and
therefore a positive theta. It will be worth slightly more as each day passes. This
is shown in Figure 7-9.
I replicated this in code here.
- the shaded area in the curve in green shows where there is negative time value for the option with dividend set to q (9% in this example)
- the blue curve is the payoff for the same option but with dividends set to 12%
- the shaded bars show the areas where Black Scholes theta is positive (blue for q = 9% and yellow for q2 = 12%) - note that the blue bar overlays parts of the yellow bar
- It closely resembles the chart from Natenberg, and given the impact of r and q on the option (and forward), this shape makes intuitive sense
- Theta itself is also closely related to that area as is visible by the coloured bars starting more or less at the intersection of the option value with the intrinsic value
It is not a formal prove but in my opinion, often intuition and charts are also very useful. Last but not least, one can also plot theta as a function of spot to see that it indeed turns positive for (deep) ITM options with positive dividends:
That also shows that Natenberg's explanation is incorrect - and it is dividends, not the depressing effect of interest rates that is at work here.
With regards to the answer provided by @nanoman, the same applies for ITM options as shown in my gifs (you need to exclude intrinsic value because it is not affected by time or IV).