Proc Timeseries and Proc Expand

Special Date and Time Interval Functions

Intck represents the number of time intervals in a given time span.
duration=intck('month',policy_sale_date,failure_date);

Intnx generates a date that is specified number of intervals from a starting value.
new_date=intnx('month',date,12);

Proc TIMESERIES
Proc timeseries analysis time-stamped transaction date with respect to trend and seasonal analysis. Time-stamped time series are often recorded at no fixed interval. Analysis often want to use time series analysis techniques that require fixed-time intervals. Therefore, the transaction data much be converted to a high frequency or accumulated to a lower frequency.

Further analysis can be performed by proc timeseries including: descriptive(global) statistics, seasonal decomposition/adjustment analysis, correlation analysis, and cross-correlation analysis.

proc timeseries data=failure_60 out=predict_series plot=(series corr decomp);
Id failure interval=month accumulate=median;
var severity;
run;

out=predict_series option specifies the resulting time series for each id.
interval=month option specifies the transaction to be accumulated on a monthly basis.
accumulate=median option specifies that the median of the transaction is to be calculated.
After the data is accumulated into a time series format, many of the procedures can be used to analyze the resulting time series data.

Proc EXPAND
proc expand is able to interpolate missing values in time series data. It provides three ways of interpolating:
1. cubic spline interpolation. It involves joining segments of third degree polynomial curves to the non-missing input values, so that the whole curve and its first and second derivatives are continuous. Method=spline in the convert statement to request this method.
2. linear spline interpolation. It fits a continuous curve by connecting successive straight line segments to the non-missing input values. Specify Method=join in the convert statement to request this method.
3. step function interpolation. It fits a discontinuous step function to the non-missing input values. Specify Method=step in the convert statement to use this method.

%let start_date=mdy(10,1,2002);
%let end_date=mdy(6,1,2011);
data structure;
length=intck('month',&start_date,&end_date);
do i=0 to length;
failure_date=intnx('month',&start_date,i);
format failure_date date.;
output;
end;
run;

proc sort data=series;
by failure_date;
run;

proc sort data=structure;
by failure_date;
run;

data series1;
merge structure(in=a) series;
by failure_date;
if a;
run;

proc expand data=series1 out=predicted_series from=month;
convert claim_amt=severity/ observed=average;
id failure_date;
run;
*used for interpolation.

convert claim_amt= lists the variables to be processed.
id option names a numeric variable that is used to identify observations in the input and output data sets.
from= option specifies the type of time interval between observations of the input data set. The interval specification can be year, semiyear, qtr, month, week, day, etc.
to= option specifies the type of time interval between observations of the output data set.
observed= option indicates the attributes of the input time series and the output series. The attribute to the following: beginning(specifies beginning of period values), middle(specifies midpoint values), end(specifies end of period values), total(specifies period totals), average(specifies period averages).
method= option specifies the method used to convert the data series. The methods supported as spline, join, step, and aggregate. 

Proc ARIMA to Fit AR(p) model

Autoregressive models are a class of models that estiamte the stochastic process underlying a time series, wher the time series values exhibit nonzero auto-correlation.

The Box-Jenkins approach to analyzing time series data is an analytical method involving three phases: identification of the process, estimation of the model parameters, and diagnostic checking.

*Identification Phase
The first phase is the identification phase. Use Proc arima by issuing the IDENTIFY statement.

proc arima data=predicted_series ;
identity var=severity scan;
run;

ACF lists the estimated auto-correlation coefficients at each lag, cov(Y_t, Y_t+i), i=1,2,.... A plot of the ACF shows the pattern of auto-correlations as the number of lags increases.

IACF is the auto-correlation function of an inverted model. the IACF of an AR is equivalent to the ACF for the same process modeled with an MA model.

PACF is the auto-correlation function at lag p accounting for the effects of all intervening observations. Thus, the PACF at lag 1 is identical to the auto-correlation at lag 1, but the PACF and ACF are different at all higher lags. Plots of the PACF reveal patterns that are useful in identifying the order of auto-regressive processes.

For an AR(p) process, the ACF plots should display the following patterns:
1. ACF tails off exponentially.
2. The PACF and IACF drop to 0 after lag p, where p is the order of the AR process. In PACF plot, if the tail-off behavior is combined with positive auto-correlations alternating with negative auto-correlation on successive lags, it is said to be contained within an exponential envelop.
3. If the ACF drops off very slowly, the time series is said to be non-stationary and will required further processing.

*Estimation Phrase
The second phrase of the ARIMA procedure is the estimation phrase.

proc arima data=predicted_series ;
identity var=severity scan;
estimate p=6;
run;

Three estimation techniques are available in Proc ARIMA:
Conditional least squares (CLS) is conditional upon the assumption that all values of the series prior to the first observation the analysis data set are equal to the mean of the series.
Unconditional least squares (ULS) is exact least squares which makes no assumption.
Maximum likelihood (ML) attempts to maximize the log of the likelihood function.

The ouput from the estimation phase lists model parameter estimates and associated significant tests, goodness of fit statistics, and other estimation diagnostics.
Use the ESTIMATE statement to invoke the estimation phrase of PROC ARIMA. An IDENTIFY statement must precede the ESTIMATE statement.

*Diagnostic Checking Phrase
The third phrase is to diagnostic checking the model.

General guidelines for checking diagnostics:
1. Check auto-correlation check of residuals. If the Q statistics are significant, you can conclude that the estimated model provides a poor fit. Q statistics is for the white noise hypothesis for the residuals of the fitted model. Each Q statistic is a chi-square produced using Ljung-Box calculation. Ljung-Box is based on the null hypothesis that the series is white noise.
2. Use the goodness-of-fit criteria(variance estimate, AIC, BIC) to compare the fit of different models estimated with same estimation method (conditional least squares, unconditional least squares, maximum likelihood).
3. Check the parameter estimates. These may indicate that model has too many parameters.
4. If the obsolete value of the last (highest-order) parameter estimates is close to 1, then the series may be non-stationary.
5. Perform a spectral analysis on the residuals to test for white noise hypothesis.

Proc arima is also used to forecast autoregressive models. ARIMA models are best for short-term forecasting because long-term forecasting is no better than using the mean of the time series.

proc arima data=predicted_series ;
identify var=severity scan;
estimate p=6;
forecast lead=60 interval=month id=failure_date out=outlead;
run;

axis1 width=1 offset=(1 pct) label=(a=90 r=0 'Severity');
axis2 width=1 offset=(1 pct) label=('Failure Date') value=(h=1.25);
symbol1 v=star ci=red height=1 cells interpol=join l=1 w=2;
symbol2 v=dot ci=green height=1 cells interpol=join l=1 w=2;
symbol3 v=dot ci=blue height=1 cells interpol=join l=2 w=1;
symbol4 v=dot ci=blue height=1 cells interpol=join l=2 w=1;
legend1 label=('') value=('Actual Claim Severity' 'Predicted Claim Severity') across=2 mode=protect position=(top center inside);
title 'Severity Prediction';
proc gplot data=outlead;
format failure_date monyy. severity dollar.;
plot (severity forecast L95 U95)*failure_date/ overlay
caxis = BLACK
ctext = BLACK
vaxis = axis1
haxis = axis2
legend= legend1
grid
hminor=0
;
run;

Proc ARIMA to Fit MA(q) model

Moving average models are a class of models that estimate the process underlying a time series, where the current time series value is related to the random errors from previous time periods. In contrast, auto-regressive models estimate a process where the current time series value is related to the actual time series values from previous time periods. Like auto-regressive processes, time series with moving average processes also exhibit nonzero auto-correlation.

The Box-Jenkins approach provides a framework for identifying the moving average process underlying a time series and for estimating model parameters. For an MA(q) process, the auto-correlation plots should display the following patterns:
1. The ACF drops to 0 after lag q.
2. The IACF and PACF tail off exponentially.

A plot of the data shows a upward trend. With this type of trend, the time series is non-stationary. A time series is said to be stationary if the mean, variance, and lag covariance are constant over time. One common way of achieving this is to difference the series. Differencing a time series involves making calculations of the type Y_t-Y_t-1 for all observations.
*take first differences to remove non-stationary trend;
data predicted_series;
set predicted_series;
severity_diff=dif(severity);
run;

The first phase is the identification phase, which you perform by issuing the IDENTIFY statement.
proc arima data=predicted_series ;
identity var=severity_diff nlag=12 scan;
run;

The second phrase of the ARIMA procedure is the estimation phrase.
proc arima data=predicted_series ;
identity var=severity_diff scan;
estimate q=2 method=ml;
run;

The ACF drops to 0 after the sixth lag, indicating MA(6).
The IACF drops to 0 after the forth lag, indicating AR(4).
The PACF has a spike at the sixth lag, indicating AR(6).

The third phrase is to forecast the time series.
proc arima data=predicted_series ;
identify var=severity_diff;
estimate q=4 method=ml;
forecast lead=60 interval=month id=failure_date out=outlead;
run;

axis1 width=1 offset=(1 pct) label=(a=90 r=0 'Severity');
axis2 width=1 offset=(1 pct) label=('Failure Date') value=(h=1.25);
symbol1 v=star ci=red height=1 cells interpol=join l=1 w=2;
symbol2 v=dot ci=green height=1 cells interpol=join l=1 w=2;
legend1 label=('') value=('Actual Claim Severity' 'Predicted Claim Severity') across=2 mode=protect position=(top center inside);
title 'Severity Prediction';
proc gplot data=outlead;
format failure_date monyy. severity dollar.;
plot (severity_diff forecast)*failure_date/ overlay
caxis = BLACK
ctext = BLACK
vaxis = axis1
haxis = axis2
legend= legend1
grid
hminor=0
;
run;
quit; 

Proc ARIMA to test Stationary.

Stationarity is a mathematical and statistical property of time series data. Much of the probability theory of time series analysis is based on the assumption that time series are stationary. If either the mean or the variance is not constant, then the time series is not stationary.

If the characteristic equation has a unit root, e.g., 1-pi1*B-pi2*B^2=1 for AR(2), it means that B=1 is a solution so then the mean drops out and the series is non-stationary. For ARIMA models, non-stationary is equivalent to one or more unit roots in the characteristic polynomial.

First, a systematic change in mean. If a series does not seem to have a constant mean (part of the definition of stationary) when graphed, that is a visible symptom of non-stationarity.

Second, a systematic change in variance. Log transformation could be able to make the variance constant. And also, log transformations are for removing exponential trends and other types of nonlinear behavior in time series data.

Third, a seasonal or periodic variation. The symptom had to do with the plot of the lag j correlation against the lag number j. Symbolizing the series value at time t by Yt, the lag j correlation is the correlation between Yt-j and Yt. If the plot of the estimated autocorrelations dies down very slowly with increasing j, this provides a second visual symptom of nonstationarity. Meanwhile, if your original time series has a sine-wave pattern, you may find it useful to use sinusoidal components as predictor variables. You can use spectural analysis to detect the presence of sinusoidal components.

There are three methods of dealing with non-stationary time series listed above: differencing, using log transformations, using sinusoidal components. Differencing is the most common method of handing non-stationary means, however, other methods are preferable for some types of series.

Differencing includes first difference, second difference, span-12 difference(for monthly data), span-4 difference(for quarterly data), etc. In many of these cases, it appeared that the series of first differences Yt - Yt-1 did seem to be stationary. Sometimes you may mistakenly take a difference when the data so not require it, which is over-differencing. If the inverse auto-correlation function (IACF) plot from the identification phase output of PROC ARIMA tails off very slowly, you may have over-differencing.

Note that if is a pure AR series, then the IACF of Yt is the same as the ACF of a pure MA process and will cut off after a few lags, i.e. the IACF behaves like the PACF of Yt .On the other hand, if is a pure MA series, then the IACF of Yt is the same as the ACF of a pure AR process and IACF will decline exponentially, i.e. the IACF again behaves like the PACF of Yt.

The Dickey-Fuller Unit Root test involves a regression of delta_Y on Y_{t-1}-Y^bar, delta_{Y-1}, ..., delta_{Y-p}, where delta indicates a first difference, i.e., Yt-Yt-1 on Y{t-1}-Y^bar, Yt-1-Yt-2, Yt-2-Yt-3,..., Yt-p-Yt-p-1. The t statistic for the parameter estimate of Y_{t-1}-Y^bar provides the statistical test of the stationary hypothesis and the distribution is call Tau. If the parameter estimate is negative and significant, the series is stationary. Otherwise, the series is not stationary, that is, it has a unit root and need to take difference.

proc arima data=predicted_series ;
i var=severity nlag=12 stationarity=(ADF=(0,1));
run;

/*
Augmented Dickey-Fuller Unit Root Tests

Type Lags Rho Pr < Rho Tau Pr < Tau F Pr > F

Zero Mean 0 0.1562 0.7175 0.11 0.7167
1 0.7651 0.8663 0.82 0.8868
Single Mean 0 -7.2238 0.2519 -2.10 0.2460 2.72 0.3774
1 -3.7934 0.5563 -1.56 0.4996 2.42 0.4546
Trend 0 -66.2436 0.0004 -6.92 <.0001 24.02 0.0010
1 -54.6745 0.0003 -5.24 0.0002 13.79 0.0010

*/

Zero Mean indicates Yt-Yt-1 on t, Yt-1, ..., with noint of Proc reg. It is assumed that the series mean is zero whether or not there are unit roots. The null hypothesis is that the time series is non-stationary, while the alternative hypothesis is that the series is stationary. Non-significant value indicates non-stationary.
Single Mean indicates Yt-Yt-1 on t, Yt-1, ..., of Proc reg with intercept. Because the model, under the alternative, has the single mean. If you simply subtracts the estimated mean for all observations, then check zero mean results.
Trend indicates Yt-Yt-1 on t, Yt-1, ..., of Proc reg with intercept.

In each case it is the coefficient of Yt-1 that indicates whether differencing needs to be done with a 0 coefficient indicating differencing. In each case, the usual t formula applies but the distribution is nonstandard. Differencing would not be appropriate as it would induce a unit root of moving average term. On the other hand, if the error is a unit root process then differencing is appropriate.

The remaining question is whether the nonstandard distribution is the same in each case. Unfortunately the answer is no, but fortunately the distributional results needed are available. For the SAS user, it is enough to know that the test computations in PROC ARIMA use the proper distributions so their p-values can be trusted. 

Proc LOGISTIC to fit Multinomial Logistic Models

Multinomial logit model or generalized logit model makes the parameters specific to the nominal. With subject-specific covariate only, the probability of response j=0..J-1 is pi_ij=exp(xi.beta)/(sum_0^(J-1))exp(xi.beta) with beta_0=0 for identification. An intercept is included in each betaj. The multinomial logit model has the property of independence from irrelevant alternatives because pi_jk=exp(xi.(beta_j-beta_k). Having too many parameters is a serious drawback of the MLM.

The multinomial logit model can be estimated in proc logistic and proc glimmix using the link=glogit.

data multinomial_response;
beta1_0 = .2;
beta1_1 = .5;
beta2_0 = -.2;
beta2_1 = .7;
seed = 1234579;
do id = 1 to 500;
x = rannor(seed);
eta1 = beta1_0 + beta1_1*x;
eta2 = beta2_0 + beta2_1*x;
prob1 = exp(eta1) / (1 + exp(eta1) + exp(eta2));
prob2 = exp(eta2) / (1 + exp(eta1) + exp(eta2));
prob3 = 1 / (1 + exp(eta1) + exp(eta2));
call rantbl(seed, prob1, prob2, prob3, y);
output;
end;
keep id x y;
run;

proc logistic data = multinomial_response;
class;
model y = x / link = glogit;
run;

Proc ARIMA to Fit ARIMA(p,d,q) model

Characteristics of Time Series Data:

1. Trend.
Trend is a long-term, consistent change in the time series values from beginning to end. It often represents a steady increase or a steady decrease in the values of time series, although other trend patterns are possible, e.g., quadratic or cubic trends in the data.

2. Seasonal Components.
Seasonal components involve a regular change in the data values that occurs at the same time in a given time period, e.g., Turkey sales are highest in November of every year.

3. Cyclical Components.
Cyclical components are long-term patterns of rising and falling data values in the time series. Many time series follow the business cycle of growth and recession in the economy.

4. Random Fluctuations.
Random fluctuations are what is 'left over' in a time series after the other three components have been accounted fro. They can represent measurement error or other unpredictable facors that affect the behaviour of a time series.

Check ACF, IACF, and PACF plots to identify p,q.
1. All three plots tail off exponentially for ARIMA(p,d,q) process.
2. For the ACF, r=max(p-1,q) beginning values followed by the behavior characteristic of an AR(p), that is, tailing off.
3. For the IACF and PACF, there are r=max(p-1,q) beginning values followed by the behavior characteristic of an MA(q), that is, tailing off.

Also, check seasonal pattern of time series as follow:
1. For seasonal AR process, the seasonal lags, s, 2s, ..., Ns, die off exponentially in ACF plot. In the IACF and PACF plots, there are spikes at each seasonal lag up to lag Ps, where there are P seasonal auto-regressive parameters necessary to model the data.
2. For seasonal MA process, the seasonal lags, s, 2s, ..., Ns, die off exponentially in IACF and PACF plot. In the ACF, there are spikes at each seasonal lag up to lag Qs, where Q is the number of seasonal moving average parameters necessary to model the data.
3. For seasonal ARMA process, the seasonal lags in the ACF tail off exponentially after lag Qs, where Q is the number of seasonal MA parameters necessary to model the data. The seasonal lags in the PACF and IACF tail of exponentially after lag Ps, where P is the number of seasonal AR parameters necessary to model the data.

proc arima data=predicted_series ;
identify var=severity(1) nlag=24 stationarity=(ADF=(0,1));
estimate p=(1,2) (12) q=(1,2) (12) method=ml;
forecast lead=60 interval=month id=failure_date out=outlead;
run;

axis1 width=1 offset=(1 pct) label=(a=90 r=0 'Severity');
axis2 width=1 offset=(1 pct) label=('Failure Date') value=(h=1.25);
symbol1 v=star ci=red height=1 cells interpol=join l=1 w=2;
symbol2 v=dot ci=green height=1 cells interpol=join l=1 w=2;
legend1 label=('') value=('Actual Claim Severity' 'Predicted Claim Severity') across=2 mode=protect position=(top center inside);
title 'Severity Prediction';
proc gplot data=outlead;
format failure_date monyy. severity dollar.;
plot (severity forecast)*failure_date/ overlay
caxis = BLACK
ctext = BLACK
vaxis = axis1
haxis = axis2
legend= legend1
grid
hminor=0
;
run;
quit; 

Proc Forecast to Automatic Forecasting.

Proc Forecast procedure is based on extrapolation of data and the forecast for a series is solely a function of the time and the past values of the series. It is a one-step method to automatically generate forecasts for hurndres of series at a time.

A time series plot is able to visually identify the presence of three possible components: level, trend, and seasonality.
First, the plot may be flat, showing only random fluctuation about some constant level.
Second, the plot may show an upward (or downward) trend over time.
Third, the plot may show some obvious seasonality.

trend= option specifies the effects of different values.:
trend=1 specifies single exponential smoothing (a constant model).
trend=2 specifies double exponential smoothing (a linear trend model).
trend=3 specifies triple exponential smoothing (a quadratic trend model).

There are three main exponential smoothing methods:
First, single exponential for level component.
Second, double exponential for trend component.
Third, stepar/winters/addwinters for seasonal component.

method=expo forecasts based on an exponential smoothing model.In fitting the trend, the parameters are allowed to change gradually over time, and earlier observations are given exponentially declining weights. It works for long-term, deterministic change.
method=stepar forecasts based on a stepwise autoregressive model. It combines time trend regression with an autoregressive model and used a stepwise method to select the lags to use for the autoregressive process. It works for short-term fluctuation.
method=winters forecasts based on a multiplicative seasonal exponential smoothing model.It combines a time trend with multiplicative seasonal factors to account for regular seasonal fluctuations in a series. It works for regular seasonal fluctuations.
method=addwinters forecasts based on an additive seasonal exponential smoothing model. It combines a time trend with additive seasonal factors to account for regular seasonal fluctuations in a series. It works for regular seasonal fluctuations.

If unspecified, the weight defaults to 1-.8^(1/trend), where trend is the value of the Trend=option.
weight=.2 for trend=1.
weight=.10557 for trend=2.
weight=.07168 for trend=3.

Single exponential smoothing and arima(0,1,1) are parallel when the moving average parameter set to smoothing constant. Because the smoothed values of St and St-1 are forecasts of Yt+1 and Yt respectively.
Double exponential smoothing and arima(0,2,2) are parallel. Proc arima constructs forecasts with the estimation suppressed and used the MA=option equals to 1-weight, 0.89442719099991 for default, in proc forecast if trend=2.
Triple exponential smoothing and arima(0,3,3) are also parallel.

proc forecast data=predicted_series interval=month method=expo trend=2 lead=70 out=predict_severity1 outfull outest=est;
id failure_date;
var severity;
run;

data severity1;
set predict_severity1;
if _type_='FORECAST';
pred1= severity;
run;

proc sort data=severity1;
by failure_date;
run;

*replicate the first observation;
data predicted_series2;
set predicted_series ;
if _n_=1 then do;
failure_date=intnx('month', failure_date,-1);
output;
failure_date=intnx('month', failure_date,1);
end;
output;
run;

*1-weight for MA;
proc arima data=predicted_series2 ;
identify var=severity(1,1) noprint;
estimate q=(1) (1) noconstant noest
ma=0.89442719099991 0.89442719099991;
forecast lead=60 interval=month id=failure_date out=predict_severity2 noprint;
run;

data severity2 ;
set predict_severity2 ;
pred2=forecast;
obs=severity;
keep failure_date pred2 obs;
run;

proc sort data=predict_severity2;
by failure_date;
run;

data mergedata;
merge severity1 severity2;
by failure_date;
run;

proc compare data=mergedata MAXPRINT=200;
var pred1;
with pred2;
id failure_date;
run;

axis1 width=1 offset=(1 pct) label=(a=90 r=0 'Severity');
axis2 width=1 offset=(1 pct) label=('Failure Date') value=(h=1.25);
symbol1 v=dot ci=red height=1 cells interpol=join l=1 w=2;
symbol2 v=dot ci=green height=1 cells interpol=join l=2 w=2;
symbol3 v=dot ci=blue height=1 cells interpol=join l=2 w=2;
legend1 label=('') across=3 mode=protect position=(top center inside);
title 'Severity Prediction';
proc gplot data=mergedata;
format failure_date monyy. ;
plot (pred1 pred2 obs)*failure_date/ overlay
caxis = BLACK
ctext = BLACK
vaxis = axis1
haxis = axis2
legend= legend1
grid
hminor=0
;
run;
quit;