@nbbo2 and @Quantuple already answered the question in their comments but if in doubt, I always think computer coding is very helpful because you can simply try it out and run a lot of calculations in one go. Since option prices and Greeks (assuming Black Scholes) are well defined, it is straightforward to set this up, and replicate the chart from the book.
Using Julia, and typing the wikipedia formulas into code, would look like this:
using Plots, Distributions, PlotThemes, DataFrames
theme(:juno)
#theme(:solarized)
Plots.scalefontsizes(1.3)
# define cdf and pdf
N(x) = cdf(Normal(0,1),x)
n(x) = pdf(Normal(0,1),x)
# generic put call pricer
function BSM(S,K,t,r,q,σ)
d1 = ( log(S/K) + (r - q + 1/2*σ^2)*t ) / (σ*sqrt(t))
d2 = d1 - σ*sqrt(t)
c = exp(-q*t)S*N(d1) - exp(-r*t)*K*N(d2)
theta = (-(S*exp(-q*t)*n(d1)*σ)/(2*sqrt(t)) + q*exp(-q*t)*S*N(d1) - r*exp(-r*t)*K*N(d2))/365
return c, theta
end
# define inputs
spot = 0.2:0.1:15 # underlying value
k = 5 # strike
r = 0.03 # interest rate: 3%
q = 0.09 # dividend = 9%
q2 = 0.12
t = 1 # time to expiry in years
σ = 0.3 # implied volatility
# compute option
opt_val = BSM.(spot,k,t,r,q,σ)
df = DataFrame(BSM.(spot,k,t,r,q,σ))
rename!(df, :2 => :theta)
theta = df[!,:theta]
df1 = DataFrame(BSM.(spot,k,t,r,q2,σ))
rename!(df1, :2 => :theta1)
theta1 = df1[!,:theta1]
x = [i < 0 ? i : 0 for i in [val[1] for val in opt_val].-[val for val in max.(spot.-k,0)]]
area = [val[1] for val in opt_val].-x
# plot option
plot(spot,[val[1] for val in opt_val],
label = "Theoretical Value with dividend = $(round(q*100, digits = 4))%",
title = "Natenberg replica",
size = (900,600),
xlabel = "Underlying Price",
ylabel = "Option Value",
linwidth = 2,
linecolor = :lightgreen)
# plot intrinsic value
plot!(spot,[val for val in max.(spot.-k,0)], label = "Intrinsic Value", legend=:topleft, linestyle = :dash)
## color the area where there is negative time value
plot!(spot, [val[1] for val in opt_val],
fillrange = area, fillalpha = 0.35, c = 3, label = "Negative Time Value with dividend = $(round(q*100,digits=4))%" )
# plot a second option with a different dividend
plot!(spot,[val[1] for val in BSM.(spot,k,t,r,q2,σ)], label = "Theoretical Value with dividend = $(round(q2*100,digits=4))%", title = "Natenberg replica")
# shade area where theta is positive
vspan!([spot[length(opt_val) - length(theta[theta .>= 0])], spot[length(spot)]],
linecolor = :grey, fillcolor = :blue, opacity = 0.05,
label ="Positive θ for dividend = $(round(q*100,digits=4))%")
vspan!([spot[length(opt_val) - length(theta1[theta1 .>= 0])], spot[length(spot)]],
linecolor = :grey, fillcolor = :yellow, opacity = 0.03,
label ="Positive θ for dividend = $(round(q2*100,digits=4))%")
It is not necessary to understand the code in full but there are a few things to note:
- the Black Scholes model and theta is the same formula as shown in Wikipedia
- 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 colored bars starting more or less at the intersection of the option value with the intrinsic value
Adding a few lines of code, similar to this answer, allows us to make this chart interactive (quality cannot be much better due to size limit):

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:

A related question for puts is here.
Edit
I wanted to comment with regards to the comments below but it turned out a bit lengthy.
I am not a 100% sure why the focus is on the depressing effect of interest rates. Copying the answer I linked above, there are 2 circumstances that can lead to the value of an european option being lower than intrinsic value
- deep ITM puts in presence of positive interest rates r>0
- deep ITM calls in presence of positive dividend yield q>0
which also coincides with the 2 circumstances under which it makes sense for an american option to be exercised early.
Insofar @nbbo2 may be right and the mixed parts of a section about puts into this section. The simplest explanation is with American puts, where stock price is zero. If the strike is 60, exercising early results in an immediate gain of 60. Due to interest rates (time value of money), receiving 60 now is preferable to later (you cannot gain more by waiting, because negative stock prices are impossible).
However, it's clearly written in terms of discounting the PV of the intrinsic value of the call, which kind of rules out that it was intended for puts. This argument is sloppy imho because it by definition rules out any time value of money. It will always be worth at least the PV of its intrinsic value (plus time value, which can be negative).
If it’s not mixing arguments for puts and calls, I think it may be because the book uses mainly intuition and simple examples as opposed to mathematical rigor. He also writes elsewhere in the book that
Theta Risk is the opposite side of gamma risk ... negative gamma
always goes hand in hand with a positive theta
Theta for deep ITM calls can be positive as we have seen. Yet, Delta will tend towards 1, which means Gamma will tend towards 0. Hence, it is perfectly reasonable to have positive Theta, with a non-negative (0) Gamma as shown below.

He also writes that delta is approximately equal to the probability that the option will finish in the money, but he moves on to explain that delta is only an approximation of the probability because interest considerations and dividends may distort this interpretation. While rates and divs will impact delta as shown here, this is not the reason delta is not the true probability. See Understanding N(d1) and N(d2): Risk-Adjusted Probabilities in the Black-Scholes Model by Lars Tyge Nielsen for a detailed explanation. Only if time to maturity and vol is low, $d1 \approx d2$ and delta will be closer to the risk-adjusted probability of the event that the option will finish in the money, which is $P(S_T > X) = N(d2)$.
While some parts may not be 100% correct, it is probably still much more useful than "correct" expositions that require a lot more detail in a book like this. It's a hands on approach with a commonsense point of view. It is sufficient to know theta can be positive, that delta is approximately equal to the probability and that usually gamma and theta have opposite signs. Anything more complicated than the concepts mentioned in the books will likely not help you much on a trading floor. You can have a look at the discussions and answers to this question to see that thinking in detail can actually be quite cumbersome, while the answer to the question likely adds no value to an option trader.
In the words of Euan Sinclair on P.XV of the preface in his book called Option Trading,
Many successful traders have had only a basic understanding of the
mathematics behind BSM and certainly have no idea of the differences
between calculus and stochastic calculus.
That's not meant to say Natenberg was "only" a successful trader but that the book was likely written in a way that builds the intuitions needed to be a successful trader, not a quant.