Question 1: Not that often, usually finite difference is use to compute it (bump and reprice). See for example this complete step by step explanation of Bloomberg's DV01 computation in SWPM. Practically, its the same though.
Question 2: No, because DV01 is an exact linear approximation to the price yield function already. That is similar to the Greek Delta in option pricing, where no one uses some sort of regression to estimate the change of the option value with respect to a change in the price of the underlying either.
TL;DR
You correctly define the price of a bond as the discounted cashflows. Insofar, I am a bit puzzled by your question. As you showed, yield and price have a direct relationship. Demand and supply changes will change price, and hence yield accordingly. A simple example with code for US T-bills can be found here.
(Modified) duration, is a mathematical derivative (rate of change) of price and measures the percentage rate of change of price with respect to yield. DV01 is the same thing measured in absolute (dollar) terms.
Let's look at the bond price formula on Wikipedia
\begin{matrix}
\left(\frac{C}{1+i}+\frac{C}{(1+i)^2}+ ... +\frac{C}{(1+i)^N}\right) + \frac{M}{(1+i)^N}
\end{matrix} \begin{matrix}
\left(\sum_{n=1}^N\frac{C}{(1+i)^n}\right) + \frac{M}{(1+i)^N}
\end{matrix}
\begin{matrix}
C\left(\frac{1-(1+i)^{-N}}{i}\right)+M(1+i)^{-N}
\end{matrix}
C = coupon payment (periodic interest payment)
N = number of payments
i = market interest rate, or required yield, or observed / appropriate yield to maturity
M = value at maturity, usually equals face value (Nominalwert)
P = market price of bond
Put in Python this looks like this
M = 1000 # nominal value
coupon = 0.05 # pct per period
c = coupon*M
i = 0.1 # YTM
N = 3 # periods
p = c*((1 - (1+i)**-N)/i) + M/(1+i)**N
yielding a price of ~875.66 in this example. Now, YTM is usally derived from a market quoted price. For this, you need some solver. It's relatively straighforward to use Newton Raphson in this example.
# Finding the yield of a coupon bearing bond
def newtonRaphson(m,p,t,c,r,f,epsilon=0.00001):
def fx(r):
cum = 0
red = m*(1+(r/100/f))**(-t*f)
for i in range(1,t*f+1):
cum += c/f*(1+(r/100/f))**(-i)
cum += red
return p-cum
def fxdx(r):
cum = 0
red = t*m*(1+(r/100/f))**-(t*f+1)
for i in range(1,t*f+1):
T = i/f
cum += t*c/f*(1+(r/100/f))**-(i+1)
cum += red
return cum
x_n = r
i = 1
while abs(fx(x_n)) > epsilon:
x_next = x_n/100 - fx(x_n)/fxdx(x_n)
x_n = x_next*100
print('Guess {}: {}'.format(i, round(x_n,5)))
i += 1
if i > 10:
print('No solution found.')
break
ytm = round(x_n,6)
return ytm
Not surprisingly, this matches the yield we used to compute the price before.
Now that we established how YTM and prices are related to each other we an look at this in a more generic way. I often find charts intuitive.
def bond_price(N,i,M,c,f):
return ((M*c/f*(1-(1+i/f)**(-f*N)))/(i/f)) + M*(1+(i/f))**(-f*N)
freq = 1
bond_price(N,i,M, coupon,freq)
freq = 2 # payments per year (semi annual)
N = 20 # maturity in years
int_changes = np.arange(start=0.001, stop=0.14, step=0.001)
M = 1000 # Nominal
c= 0.05 # 5% coupon
[plt.plot(int_changes, bond_price(i ,int_changes,M,c, freq), label = f'Bond with {i} year(s) maturity') for i in np.arange(start=1, stop=32, step=10)]
plt.xlabel("YTM")
plt.ylabel("Bond price")
plt.title('Bond prices as a function of interest rates')
plt.axvline(x=0.05, color='black', linestyle='-.', linewidth = 0.4)
plt.axvline(x=0.02, color='black', linestyle='-.', linewidth = 0.4)
plt.axhline(y=1688, color='c', linestyle='-', linewidth = 0.4)
plt.legend(prop={'size': 10})
plt.show()
Duration is just a linear approximation to this curve, similar to how delta works for options (see here). Implementing the formula in Wikipedia you compute it as such.
M = 100
c = 0.07
freq = 1
N = 5
r = 0.07
year = [i for i in range(1,N+1)]
Cashflow = [round(M * c, 3) if i < N else 100 + round(M * c, 3) for i in range(1,N+1)]
dcf_t = [f"1/((1+{r})^{i})" for i in range(1,N+1)]
dcf = [round(1/((1+r)**i), 5) for i in range(1,N+1)]
cf_weighted = [i*round(M * c, 3) * round(1/((1+r)**i), 5) if i < N else i*(100 + round(M * c, 3)) * round(1/((1+r)**i), 5)for i in range(1,N+1) ]
cf_disc = [round(M * c, 3) * round(1/((1+r)**i), 5) if i < N else (100 + round(M * c, 3)) * round(1/((1+r)**i), 5)for i in range(1,N+1) ]
df = pd.DataFrame({"Year" : year, "Cashflow" : Cashflow, "Discount Formula" : dcf_t, "Discount Factor" : dcf, "NPV Weighted CF" : cf_weighted, "NPV CF" : cf_disc})
df
Using a free online bond calculator or WolframAlpha provides identical values.
Computing MD via finite difference (similar to how BBG computes DV01 and the delta example linked above) works like this.
$$ Modified Duration = \frac{P^{-} - P^{+}}{2*P_0*dY}$$
where
$dY$ is difference in the interest rate / YTM (frequently 100bp)
$P^{-}$ Price of bond if rates down
$P^{+}$ Price of bond if rates up
$P_0$ Market value of bond
Adding duration to the chart above looks like this
int_changes = np.arange(start=0.001, stop=0.14, step=0.001)
plt.plot(int_changes, bond_price(N ,int_changes,M,c, freq), label = f'Bond with {N} year(s) maturity')
plt.plot(int_changes, [bond_price(N,r,M, c,freq) - duration/(1+r)*i*100 + duration/(1+r)*r*100 for i in int_changes], label = f'Duration approximation', linestyle='dashdot')
plt.xlabel("YTM")
plt.ylabel("Bond price")
plt.title('Bond prices as a function of interest rates')
plt.axvline(x=r, color='black', linestyle='-.', linewidth = 0.4)
plt.axhline(y=bond_price(N,r,M, c,freq), color='c', linestyle='-', linewidth = 0.4, label = f'Market price for YTM = {r}' )
plt.legend(prop={'size': 10})
plt.show()
Question 2:
So what about your linear regression idea? You are essentially proposing to use a linear functional form for a convex function. Our example already computed a few observation. Therefore, we can quickly use statsmodel to write the following Python code which implements your idea with our example from above. As you can see, the slope is actually quite similar to the duration we just computed (in absolute values).
x = int_changes
y = bond_price(N ,int_changes,M,c, freq)
model = sm.OLS(y, sm.add_constant(x))
results = model.fit()
results.summary()
Charting this with a scatter plot also shows the similarity (just add plt.axline(xy1=(0, b), slope=m, label=f'$y = {b:+.1f} {m:.1f}x $', color='r')
to the chart.
The more pronounced convexity (e.g. the longer maturity), the worse your idea gets.
As you can see, for a 30 year bond, your linear regression logic is misstating the change in the current price of the bond be quite a bit, whereas DVO1 (duration) is still an exact linear approximation at the point of interest (the current market price).
Of course, you can change the functional form and make the OLS regression fit the curve, but that will really just provide a very similar result to using duration, and adding convexity (or using Delta and adding Gamma in the options example I linked), albeit in a more convoluted way, and adding estimation errors.
