Investigating variability at the earliest stages of low-mass star formation is fundamental in understanding how a protostar assembles mass. While many simulations of protostellar disks predict non-steady accretion onto protostars, deeper investigation requires robust observational constraints on the frequency and amplitude of variability events characterized across the observable SED. In this study, we develop methods to robustly analyze repeated observations of an area of the sky for submillimeter variability in order to determine constraints on the magnitude and frequency of deeply embedded protostars. We compare 850 m JCMT Transient Survey data with archival JCMT Gould Belt Survey data to investigate variability over 2-4 year timescales. Out of 175 bright, independent emission sources identified in the overlapping fields, we find seven variable candidates, five of which we classify as Strong, and the remaining two we classify as Extended to indicate that the latter are associated with larger-scale structure. For the Strong variable candidates, we find an average fractional peak brightness change per year of 4.0 % yr − 1 , with a standard deviation of 2.7 % yr − 1 . In total, 7% of the protostars associated with 850 m emission in our sample show signs of variability. Four of the five Strong sources are associated with a known protostar. The remaining source is a good follow-up target for an object that is anticipated to contain an enshrouded, deeply embedded protostar. In addition, we estimate the 850 m periodicity of the submillimeter variable source, EC 53, to be 567 32 days, based on the archival Gould Belt Survey data.