The advent of horizontal drilling technology and hydraulic fracture completions has fundamentally altered the geometry of flow in shale gas systems, requiring new decline curves to forecast production. In this work we perform automatic decline analysis on Marcellus gas shale wells, and predict ultimate recovery for each well. We use a minimal model that captures the basic physics and geometry of the extraction process. This model allows us to scale all individual wells' gas productions to fit upon a single, universal curve. This scaling curve predicts the production rate will decline as 1 divided by the square root of time early in the life of the well, and eventually rate declines exponentially. We fit cumulative production to the results of this model to determine time to boundary-dominated flow and ultimate recovery. For wells that have not entered boundary dominated flow, we use a new technique to estimate ultimate recovery, relying on early production. We find 404 wells that are experiencing boundary-dominated flow. One would expect the time to boundary dominated flow varies due to well completion and reservoir properties. Instead, we find ultimate recovery in wells that have experienced boundary dominated flow is largely predicted by initial production, and time to exponential decline varies surprisingly little. Utilizing this information, we forecast gas production for 5,275 horizontal wells in Pennsylvania and West Virginia for a period of 25 years. We find that the median 25-year EUR is 3.6 Bcf with a standard deviation of 3.6 Bcf. The median time to interference is 3.9 years. From these results, we propose a new method for estimating time to boundary dominated flow for young wells. Our key discovery is that wells can have their EURs estimated early in life with surprising accuracy.