Precession in an accretion-powered pulsar is expected to produce characteristic variations in the pulse properties. Assuming surface intensity maps with one and two hotspots, we compute theoretically the periodic modulation of the mean flux, pulse-phase residuals and fractional amplitudes of the first and second harmonic of the pulse profiles. These quantities are characterised in terms of their relative precession phase offsets. We then search for these signatures in 37 days of X-ray timing data from the accreting millisecond pulsar XTE J1814-338. We analyse a 12.2-d modulation observed previously and show that it is consistent with a freely precessing neutron star only if the inclination angle is < 0.1 degrees, an a priori unlikely orientation. We conclude that if the observed flux variations are due to precession, our model incompletely describes the relative precession phase offsets (e.g. the surface intensity map is over-simplified). We are still able to place an upper limit on ε of 3.0 x 10^{-9} independently of our model, and estimate the phase-independent tilt angle \theta; to lie roughly between 5 and 10 degrees. On the other hand, if the observed flux variations are not due to precession, the detected signal serves as a firm upper limit for any underlying precession signal. We then place an upper limit on the product ε cos(\theta) of ≤ 9.9 x 10^{-10}. The first scenario translates into a maximum gravitational wave strain of 10^{-27} from XTE J1814-338 (assuming a distance of 8 kpc), and a corresponding signal-to-noise ratio of ≤ 10^{-3} (for a 120 day integration time) for the advanced LIGO ground-based gravitational wave detector.