We investigate crustal deformation patterns associated with magmatic activity at Taal Volcano, an active tholeiitic volcano located within the Taal caldera in SW Luzon, Philippines. We compare elastic and combined elastic-viscoelastic modeling results with deformation observed by GPS on and near the volcano. Continuous dual-frequency and single-frequency (L1) GPS data between 1998 to 2005 provide evidence for smoothly transitioning periods of inflation and deflation centered under the volcano. These deformation phases typically last 3-9 months, with rates exceeding 200 mm/yr. The most prominent observed inflationary phases took place in February-November 2000 and June 2004 to March 2005. The 2000 period of inflation was characterized by approximately 120 mm of uplift of the center of Volcano Island relative to the northern caldera rim. The 2004-05 inflation, meanwhile, was characterized by ~73 mm/yr extension of the volcanic edifice and ~50 mm uplift. Observed surface deformation is used for determining the best fitting Mogi point source, presumably associated with magmatic intrusion. Elastic modeling of the source region indicates a relatively stable point source, located under Volcano Island's central crater, at depths of ~5 km with time varying pressure or volume changes. We investigate finite element models of volcanic deformation using axisymmetric geometry, including a centrally located small spherical source within an elastic half space. Models were tested using annuli of viscoelastic shells, embedded in multi-layered elastic half-space. Using simply varying pressure histories as input, a variety of forward models are fit to the time history of continuously observed GPS deformation for two stations on volcano island relative to those situated outside the caldera. Though the precise geometry, volume or material properties that surround the presumably magmatic pressure source are non-unique, they use more geologically appropriate material properties (including viscoelasticity), which allow for a better understanding of the magmatic source processes. The use of viscoelastic models allow for simpler pressure histories and require significantly reduced overall pressure increases, relative to equivalently fitting purely elastic models.