Total Harmonic Distortion calculation and its origins

I ve been researching about harmonic distortion and THD. From what i ve seen on the internet, HD is caused by non-linear signal processing. I want to calculate it using an digital signal processor. For this, i suppose i need to compare the distorted signal with the original one.

My question is, is it possible to generate THD in a digital signal only by digital processing? Or i need to cause it in an analogic enviroment then convert to digital for detection?

Cause i ve only found examples of it being caused by analogic processing like amplifiers.

y(t) = a0 + a1*x(t) + a2*x(t).*x(t) + a3*x(t).*x(t).*x(t) + ....

where y(t) is the system output, x(t) is the input, and an are coefficients that depend on the system you're modeling.
As an example, you can do a simple simulation of mild harmonic distortion by setting x(t) equal to the a cosine signal, and setting a1=2, a3=-0.1, and all other coefficients equal to 0. a1 is the system gain, and a3 is the non-linear gain. By setting a3 to larger negative values, you increase the nonlinear effect -- you can think of it as an amplifier getting closer to saturation.
If you plot the spectra of x(t) and y(t), you'll clearly see harmonic distortion occuring.