The Hilbert is very simple to implement with high accuracy for narrow band signals such as the OP presents (we just need to introduce a perfect quadrature phase shift and maintain the amplitude balance-and with a single tone as the OP has done, this can be accomplished perfectly down to our numerical precision). But for wideband signals the Hilbert can be more challenging depending on accuracy required.
A dominant reason for time ripple is delay mismatch (in particular in the OP’s case there would be an additional half sample delay introduced in the assumed delay of the Hilbert filter). This would occur in this case if the length $N$ of the filter is 64 since the length is one higher than the order, and the delay is $(N-1)/2$. Filter truncation will cause ripple in frequency, however not in time- ripples in time will be introduced by quadrature and amplitude mismatch, and the delay leads to quadrature or phase error. For this simple case of a narrow band tone in particular, amplitude and phase matching can be accomplished with even shorter length filters.
The Hilbert as computed is half the length of the filter, and since an odd length filter is used, there will be a half sample residual delay. For cases when we can sample at twice the rate, the proper half sample delay correction can be introduced by delaying by an integer number of samples at the higher rate and then after decimating by two will result will be a perfect half sample delay at all frequencies.
To understand how the half sample delay results in the ripple (and similarly if the amplitude was not balanced), note that the half sample delay will be a phase shift at any given frequency, and from the equations and graphics below we can see how ripple in the time domain would be introduced:
First consider the analytic signal (which is the sum of the signal plus the Hilbert transform):
$$x_a(t) = x(t) + j\hat x(t)$$
For a sinusoid this is given as:
$$x_a(t) = \cos(\omega t) + j\sin (\omega t)$$
And from this case with Euler's equation and the graphic below, we see how this reduces to $x_a(t) = e^{j\omega t}$; a single rotating phasor on the complex plane whose magnitude (envelope) is always 1.
Now if we were to delay $x(t)$ relative to $\hat x(t)$, this would be a fixed phase at one particular frequency, and the result would be two phasors in counter rotation, and for small angles it would appear as in the graphic below with a larger phasor added to a smaller phasor in counter rotation (consistent with a half sample offset of a 1 KHz sinusoid at a sampling rate of 16 KHz, consider the normalized radian frequency is $\omega = 2\pi(1e3)/16e3 = 0.392$ radians per sample, so a half sample delay in the OP's case is 0.196 radians).
The formula in the graphic is derived as follows using $y(t)$ as our faulty analytic signal with phase introduced:
$$y(t) = \cos(\omega t + \phi) + j \sin (\omega t + \phi)$$
Expanding into exponentials using Euler's formula:
$$= \frac{e^{j(\omega t + \phi)}+e^{-j(\omega t + \phi)}}{2} + j\frac{e^{j\omega t}-e^{-j\omega t}}{2j}$$
$$ = \frac{e^{j\phi/2}e^{j(\omega t + \phi/2)}+e^{-j\phi/2}e^{-j(\omega t + \phi/2)}}{2} + \frac{e^{-j\phi/2}e^{j(\omega t + \phi/2)}-e^{j\phi/2}e^{-j(\omega t+\phi/2)}}{2}$$
$$ = \frac{e^{j\phi/2}+e^{-j\phi/2}}{2} e^{j(\omega t + \phi/2)} + \frac{e^{-j\phi/2}-e^{j\phi/2}}{2} e^{-j(\omega t + \phi/2)}$$
$$ = \cos(\phi/2)e^{j(\omega t + \phi/2)}-j\sin(\phi/2)e^{-j(\omega t + \phi/2)}$$
If we were to observe the smaller phasor relative to the larger phasor we immediately see how the addition of the two phasors results in envelope ripple in time! This can be compared to the actual results to see if delay is the only factor (or if there was an amplitude error as well).
Note: An FFT on the data would reveal that this is the case: the analytic signal has positive frequency components only (meaning all bins above N/2 should be zero). Quadrature imbalance in phase and amplitude causes "side-band leakage" which means the negative frequencies (those bins above N/2) will be non-zero. This is indeed how the Matlab hilbert function works (which returns the analytic signal, not actually the Hilbert transform alone but the imaginary component will be the Hilbert): the hilbert function under the hood takes an FFT of the data and zeros the negative frequencies and returns the inverse FFT of that.
