Temperature coefficient of diode
Preceding the question

As I'm sure you can find repeated in many places, the simple Shockley diode equation is:

$$I_\text{D}=I_{\text{SAT}\left(T\right)}\left(e^\frac{V_\text{D}}{\eta \:V_T}-1\right), \text{where }V_T=\frac{k\: T}{q}$$

The emission coefficient, \$\eta\$, is usually taken to be 1, by default. The saturation current \$I_{\text{SAT}\left(T\right)}\$ is effectively a \$y\$-axis intercept point found by projecting the log-chart plot of diode current vs diode voltage. The -1 term in the equation removes this offset and brings the results into expectation where the diode voltage is \$0\:\text{V}\$.

Solving for \$V_\text{D}\$:

$$V_\text{D}=\eta\:V_T\:\operatorname{ln}\left(\frac{I_\text{D}} {I_{\text{SAT}\left(T\right)}}+1\right)\approx \frac{k\: T}{q}\:\operatorname{ln}\left(\frac{I_\text{D}} {I_{\text{SAT}\left(T\right)}}\right)$$

Just as you write.

Your question is about why, if this is all there is, that there the sign of the derivative of diode voltage \$V_\text{D}\$ with respect to temperature \$T\$ is negative, rather than positive.

It's a good question.

Answer

The missing element is that the commonly found Shockley equation is short-hand. (I've been hinting about it, above, with the use of \$I_{\text{SAT}\left(T\right)}\$, earlier.)

As it turns out, the temperature dependence of \$I_{\text{SAT}\left(T\right)}\$ is huge. The approximate equation looks like this:

$$I_{\text{SAT}\left(T\right)}=I_{\text{SAT}\left(T_\text{nom}\right)}\cdot\left[\left(\frac{T}{T_\text{nom}}\right)^{3}e^{\frac{E_g}{k}\cdot\left(\frac{1}{T_\text{nom}}-\frac{1}{T}\right)}\right]$$

\$E_g\$ is the effective energy gap (in eV) and \$k\$ is Boltzmann's constant (in appropriate units.) \$T_\text{nom}\$ is the temperature at which the equation was calibrated, of course, and \$I_{\text{SAT}\left(T_\text{nom}\right)}\$ is the extrapolated saturation current at that calibration temperature.

This formula heavily depends upon fundamental thermodynamics theory and the Boltzmann factor, which you can easily look up and is above represented by the factor: \$e^{\frac{E_g}{k}\cdot\left(\frac{1}{T_\text{nom}}-\frac{1}{T}\right)}\$. It's based on the simple ratio of the numbers of states at different temperatures; really no more complex than fair dice used in elementary probability theory. Perhaps the best introduction to the Boltzmann factor is C. Kittel, "Thermal Physics", John Wiley &amp; Sons, 1969, chapters 1-6 in particular.

The power of 3 used in the equation above is actually a problem, because of the temperature dependence of diffusivity, \$\frac{k T}{q} \mu_T\$. And even that, itself, ignores the bandgap narrowing caused by heavy doping. In practice, the power of 3 is itself turned into a model parameter.

The power of 3 is approximate and is sometimes treated as an "adjustable" parameter. But close.

Summary

Once you realize the fuller form of the Shockley equation, you are now able to work out the details about why the temperature variation is so different -- not only in magnitude, but also in sign. The thermal voltage component would normally imply a positive sign. But the \$I_{\text{SAT}\left(T\right)}\$ factor overwhelms this response and changes the sign, itself.