The coupled perturbed Kohn-Sham (CPKS) computational scheme for the evaluation of electric susceptibility tensors in periodic systems, recently implemented in the CRYSTAL code, has been extended to third-order. It is, then, used to obtain static electronic hyperpolarizabilities of zigzag BN nanotubes for the first time. This procedure, which is fully analytic in all key steps, requires a double self-consistent treatment for taking into account the first- and second-order response of the system to the applied field. The performance of different functionals is compared and the B3LYP hybrid is ultimately chosen for calculations on nanotubes having radii as large as R = 20 angstrom (6-200 atoms in the unit cell). Such large radii are sufficient to give the pure longitudinal component of the (hyper)polarizability tensors to within 1% of the "exact" hexagonal BN monolayer limit. Other tensor components involving the transverse direction converge more slowly. They can, however, be extrapolated to the monolayer limit to within 4% accuracy except for the pure transverse second hyperpolarizability, which has an error of 13% in that limit.