How to simulate a Lorentzian broadened spectrum in MATLAB?
Update : In the given data by the OP, the x scale is not evenly spaced so the discrete convolution approach will not work to generate a spectrum until and unless your data is evenly spaced! The trick mentioned by EdV is a very nice shortcut. Here, if you really want to play around with convolution via MATLAB, then read below.
Solution: You can satisfy yourself as follows:

Make sure you have a even spacing of the x-axis. Your uneven axis threw me off, and it took a while to discover why it is not working. I wish you mentioned it earlier.
Your sticks are basically scaled delta functions as Prof. Ed said.
Generate a vector for a single Lorentzian peak, whose center is zero, which is also evenly spaced. The x-axis of the Lorentzian would be symmetric [-length (wavenumber/2):Sampling rate:length(wavenumber/2)-Sampling rate].
Perform convolution of your spectrum and this Lorentzian.

Let us start with a Dirac delta centered at point $x.$ I don't know how to draw a displace Dirac delta in MATLAB. I manually made it in an Excel file. If you convolute it with a zero centered Lorentzian, then you will get a Lorentzian centered at $x.$ This is called the sifting property.
Basically, all you have to ensure is that your Lorentzian is centered at zero. First, play with convolution using single stick. Once you correct the code, so that the line position does not change, apply it on the entire spectrum.

A great book is Bracewell's Fourier Transform And Its Applications. It is all there. From MIT 2.14 / 2.140 Analysis and Design of Feedback Control Systems, Spring 2007 handout:

1.1 The “Sifting” Property of the Impulse
When an impulse appears in a product within an integrand, it has the property of ”sifting” out the value of the integrand at the point of its occurrence:
$$\int_{-\infty}^\infty f(t)δ(t - a)\,\mathrm dt = f(a)\tag{5}$$
This is easily seen by noting that $δ(t − a)$ is zero except at $t = a,$ and for its infinitesimal duration $f(t)$ may be considered a constant and taken outside the integral, so that
$$\int_{-\infty}^\infty f(t)δ(t - a)\,\mathrm dt = f(a)\int_{-\infty}^\infty δ(t - a)\,\mathrm dt = f(a)\tag{6}$$
from the unit area property.