The existence of a large population of Compton thick (CT, $N_{H}>10^{24} cm^{-2}$) AGN is a key ingredient of most Cosmic X-ray background synthesis models. However, direct identification of these sources, especially at high redshift, is difficult due to the flux suppression and complex spectral shape produced by CT obscuration. We explored the Chandra COSMOS Legacy point source catalog, comprising 1855 sources, to select via X-ray spectroscopy, a large sample of CT candidates at high redshift. Adopting a physical model to reproduce the toroidal absorber, and a Monte-Carlo sampling method, we selected 67 individual sources with >5% probability of being CT, in the redshift range $0.04<z<3.5$. The sum of the probabilities above $N_{H}>10^{24} cm^{-2}$, gives a total of 41.9 effective CT, corrected for classification bias. We derive number counts in the 2-10 keV band in three redshift bins. The observed logN-logS is consistent with an increase of the intrinsic CT fraction ($f_{CT}$) from $sim0.30$ to $sim0.55$ from low to high redshift. When rescaled to a common luminosity (log(L$_{rm X}$/erg/s)$=44.5$) we find an increase from $f_{CT}=0.19_{-0.06}^{+0.07}$ to $f_{CT}=0.30_{-0.08}^{+0.10}$ and $f_{CT}=0.49_{-0.11}^{+0.12}$ from low to high z. This evolution can be parametrized as $f_{CT}=0.11_{-0.04}^{+0.05}(1+z)^{1.11pm0.13}$. Thanks to HST-ACS deep imaging, we find that the fraction of CT AGN in mergers/interacting systems increases with luminosity and redshift and is significantly higher than for non-CT AGN hosts.