College Enrollment Forecast (State Level)
Using ARIMA to forecast college enrollment for 2015 and 2016 at state level.
Postsecondary School Enrollment Forecast (State Level)¶
NCES publishes its school enrollment data every year here: https://nces.ed.gov/ipeds/datacenter/
However, the latest you data can find is always 1 year behind, in order to get the latest estimate, below shows my implementation using ARIMA to forecast current and next year's enrollment number.
NCES has forecast on national level and few demographic attribute, however, there's no state level forecast.
1. Aggregate College Enrollment by State¶
In [3]:
# load package
import matplotlib.pyplot as plt
%matplotlib inline
import pandas as pd
from statsmodels import api as sm
import numpy as np
# turn off warning
import warnings
warnings.filterwarnings('ignore')
# load data
enrollment = pd.read_excel('DATA_SOURCE\ENROLLMENT_FORECAST_DS.xlsx', sheetname='EnrollmentByInstitution')
inst_state = pd.read_excel('DATA_SOURCE\ENROLLMENT_FORECAST_DS.xlsx', sheetname='Institution')
inst_state.columns = ['UnitId', 'Name', 'City', 'State']
state_enrollment = enrollment.merge(inst_state.drop('Name', axis = 1))
# a preview
state_enrollment.head(2)
Out[3]:
Transpose the dataframe, use state as column and convert to dateindex
In [5]:
# some plastic surgery on the df
# group sum by state
to_group = state_enrollment.drop(['UnitId', 'Institution Name', 'City'], axis = 1)
to_group = to_group.replace('-', 0)
state_enrollment_grp = to_group.groupby('State')
df = state_enrollment_grp.agg(sum).reset_index()
df = df.transpose()
# columns as state
df.columns = df.ix[0,:].values
df = df.drop(['State'])
# add time series index
df.index = [pd.Period(x, freq = 'A') for x in range(1990,2015)]
df.head()
Out[5]:
2. Prepare the Model¶
Timeseries Model: ARIMA - Details on ARIMA
Error Metrics: MAPE - Details on MAPE
In [6]:
def produce_forecast(hold_out_n, inog, order, exog = None):
"""
CV module
input:
hold_out_n (hold out n years from training for cross validation)
inog: endogenous variable (the enrollment number to forecast)
exog: exogenous variable (other variable to forecast, if any)
order: (p, d, q) parameter of ARIMA model
output:
print (average MAPE across all states)
mape by state
"""
def forecast_state(inog, exog, order, hold_out_n):
# prep exog
if exog is not None:
train_exog = exog[:-hold_out_n]
hold_out_exog = exog[-hold_out_n:]
else:
train_exog = hold_out_exog = None
# forecast wrapper
m = sm.tsa.ARIMA(inog, order=order, exog=train_exog)
model = m.fit()
return model.forecast(hold_out_n, exog = hold_out_exog, alpha = .95)[0]
forecast_y = []
# data needed
Y = inog.ix[:-hold_out_n,]
yTrue = inog.ix[-hold_out_n:,]
if exog is not None:
X = exog
else:
X = None
for state in Y.columns:
stateY = Y[state]
yHat = forecast_state(stateY.values.astype('float64'), X, order, hold_out_n)
forecast_y.append(yHat)
test = yTrue
yHat = np.array(forecast_y).transpose()
error_matrix = abs(yTrue-yHat)/yHat
print('MAPE for ARIMA %i, %i, %i' % order)
print((np.sum(error_matrix, axis = 1)/Y.shape[1]))
plt.figure(figsize = (15,3))
plt.title('mape');
plt.plot(error_matrix.ix[0], '-o');
plt.plot(error_matrix.ix[1], '-o');
plt.xticks(range(Y.shape[1]), df.columns, rotation = 'vertical');
plt.legend(error_matrix.index)
return error_matrix
In [7]:
error_100 = produce_forecast(2, df, (1,0,0), exog = None)
Looks good!
3. Add National Umemployment Rate as Exog¶
In [8]:
# load umemployment data
umemployment = pd.read_excel('DATA_SOURCE\ENROLLMENT_FORECAST_DS.xlsx', sheetname='Unemployment Rate')
umemployment.index = [pd.Period(x, freq = 'A') for x in umemployment.Year]
umemployment.drop(['Year'], axis = 1, inplace=True)
In [11]:
# umemployment lag 1
error_100_ue1 = produce_forecast(2, df, (1,0,0), exog = umemployment.ix[33:-1,0])
In [13]:
# umemployment lag 0 + lag 1
error_100_ue01 = produce_forecast(2, df, (1,0,0), exog = np.vstack([umemployment.ix[33:-1,0], umemployment.ix[34:,0]]).transpose())
Slight improvement.
4. Add Disposable Income as Exog¶
In [14]:
# load disposable income data
dp = pd.read_excel('DATA_SOURCE\ENROLLMENT_FORECAST_DS.xlsx', sheetname='DISPOSABLE_INCOME')
dp.columns = ['Year', 'DSPI']
dpg = dp.groupby('Year')
dp = dpg.mean().reset_index()
dp.index = [pd.Period(x, freq = 'A') for x in dp.Year]
dp.drop(['Year'], axis = 1, inplace=True)
In [17]:
# umemployment lag 0 + lag 1 + lag 5 disposable income
# not helping, but lag 5 has the lowest avg mape
exog = np.vstack([umemployment.ix[33:-1,0], umemployment.ix[34:,0], dp.ix[27:-5,0]]).transpose()
error_100_ue01 = produce_forecast(2, df, (1,0,0), exog = exog)
Write Forecast File¶
In [18]:
# 2015/6's umpmployment rate forecast from fed = [5.0, 4.8]
# http://www.federalreserve.gov/monetarypolicy/files/fomcprojtabl20150917.pdf
umemployment.ix[pd.Period(2015, freq = 'A')] = 5.0
umemployment.ix[pd.Period(2016, freq = 'A')] = 4.8
exog = np.vstack([umemployment.ix[34:,0], umemployment.ix[33:-1,0]]).transpose()
order = (1,0,0)
forecast_y = []
for state in df.columns:
stateY = df[state]
# yHat = forecast_state(stateY.values.astype('float64'), X, order, hold_out_n)
m = sm.tsa.ARIMA(stateY.values.astype('float64'), order=order, exog=exog[:-2])
model = m.fit()
yHat = model.forecast(2, exog = exog[-2:], alpha = .95)[0]
forecast_y.append(yHat)
In [19]:
output = pd.DataFrame(forecast_y)
output = output.transpose()
output.columns = df.columns
output.index = [pd.Period(2015, freq = 'A'), pd.Period(2016, freq = 'A')]
pd.concat([df, output]).transpose().head()
# pd.concat([df, output]).transpose().to_csv('output.csv')
Out[19]: