A systematic {it ab initio} study of the second-order nonlinear optical properties of BN nanotubes within density functional theory in the local density approximation has been performed. Highly accurate full-potential projector augmented-wave method was used. Specifically, the second-harmonic generation ($chi_{abc}^{(2)}$) and linear electro-optical ($r_{abc}$) coefficients of a large number of the single-walled zigzag, armchair and chiral BN nanotubes (BN-NT) as well as the double-walled zigzag (12,0)@(20,0) BN nanotube and the single-walled zigzag (12,0) BN-NT bundle have been calculated. Importantly, unlike carbon nanotubes, both the zigzag and chiral BN-NTs are found to exhibit large second-order nonlinear optical behavior with the $chi_{abc}^{(2)}$ and $r_{abc}$ coefficients being up to thirty times larger than that of bulk BN in both zinc-blende and wurtzite structures, indicating that BN-NTs are promising materials for nonlinear optical and opto-electric applications. Though the interwall interaction in the double-walled BN-NTs is found to reduce the second-order nonlinear optical coefficients significantly, the interwall interaction in the single-walled BN-NT bundle has essentially no effect on the nonlinear optical properties. The prominant features in the spectra of $chi_{abc}^{(2)}(-2omega,omega,omega)$ of the BN-NTs are suscessfully correlated with the features in the linear optical dielectric function $epsilon (omega)$ in terms of single-photon and two-photon resonances.