Bose-Hubbard model is mainly used to describe and study the interactions between neutral atomic gases trapped in an optical lattice [1] and Josephson junction arrays [2]. It is one of the toy models to understand quantum phase transitions, i.e. a phase transition exists between the Mott insulator state and the super- fluid state. Analytical solutions are limited to obtaining the ground state energy for small systems, whereas, computational studies can be done for larger system sizes. We applied the variational Monte Carlo method to the Bose-Hubbard model based on projected wavefunctions, i.e. Baeriswyl and Gutzwiller-Baeriswyl. Even though our method can be applicable to any dimension, we only consider the one dimensional case in this thesis. We expressed observables in forms of averages over configurations to which we can apply Monte Carlo sampling techniques. Our results for both Baeriswyl and Gutzwiller projections are in qualitatively good agreement with the known calculations of the phase diagram [3,4]. Furthermore, we introduced a new method, apart from other known methods [5, 6], based on the Drude weight [7–9] to calculate the superfluid fraction, which can also be extended to observe BCS superconductivity [10].