We present a method to investigate the kinetics of protein folding and the dynamics underlying the formation of secondary and tertiary structures during the entire reaction. By writing the solution of the Fokker-Planck equation in terms of a path integral, we derive a Hamilton-Jacobi variational principle from which we are able to compute the most probable pathway of folding. The method is applied to the folding of the Villin headpiece subdomain simulated using a Go model. An initial collapsing phase driven by the initial configuration is followed by a rearrangement phase, in which secondary structures are formed and all computed paths display strong similarities. This completely general method does not require the prior knowledge of any reaction coordinate and is an efficient tool to perform simulations of the entire folding process with available computers.