An efficient continuous-time path-integral Quantum Monte Carlo algorithm for the lattice polaron is presented. It is based on Feynman's integration of phonons and subsequent simulation of the resulting single-particle self-interacting system. The method is free from the finite-size and finitetime-step errors and works in any dimensionality and for any range of electron-phonon interaction. The ground-state energy and effective mass of the polaron are calculated for several models. The polaron spectrum can be measured directly by Monte Carlo, which is of general interest.PACS numbers: 71.38.+i, 02.70.LqThe past few years have witnessed a rapid development of continuous-time Quantum Monte Carlo (QMC) algorithms for quantum-mechanical lattice models. The driving force behind this is the desire to eliminate the systematic errors introduced by the Trotter decomposition in the standard discrete-time QMC methods [1]. The main idea then is to regard the imaginary-time evolution of a particle or spin configuration as a continuoustime Poisson process with "events" being either a particle jump or a spin flip. In this way continuous-time QMC algorithms were developed for a particle in an external potential [2], Heisenberg model [3,4], t − J model [5], bosonic Hubbard model [6], and Fröhlich polaron [7].In this Letter I present a continuous-time path-integral QMC algorithm for the lattice polaron, i.e., an electron strongly interacting with phonons on a lattice. The method combines analytical integration of the phonon degrees of freedom with the advantages of the continuoustime formulation of the Monte Carlo process. The method is universal. It works for infinite lattices in any dimensionality and for any radius of electron-phonon interaction. It is also free from the systematic finite-timestep errors which were an undesirable feature of the original (discrete-time) QMC algorithm based on the integration of phonons [8]. In fact, the new method allows for exact (in the QMC sense) calculation of the ground-state energy, effective mass, and even spectrum of the polaron. The numerical accuracy of the method is 0.1−0.3%, which is not as good as that of exact diagonalization [9,10] and density-matrix renormalization group [11] schemes, but is good enough for practical purposes.I begin with the important question of which quantities can be calculated with a path-integral QMC method. The traditional scheme samples paths which are periodic in imaginary time, see, e.g., Refs. [1,8]. This allows for the calculation of thermodynamic properties such as internal energy or specific heat, as well as some static correlation functions, but not dynamic properties of the system. In Ref. [12] it was shown that the polaron effective mass, an important dynamic characteristic, can be measured on the ensemble of paths with open boundary conditions in imaginary time (BCIT). Here I extend the argument to the whole polaron spectrum. Consider two quantities. The first one, Z P = i i|e −βH |i δ Pi,P , is by definition the partition function of many-body s...