This is best done through simulation. See my MATLAB code example and explanation below:
%% Get S&P 500 price series
d=fetch(yahoo,'^GSPC','Adj Close','1-jan-2014','30-dec-2014');
n = 1; % # of shares
p = d(end:-1:1,2); % share price, the dates are backwards
PV0 = n*p(end); % portfolio value today
%%
r=price2ret(p,[],'Continuous'); % get the continous compounding returns
Mdl = arima('ARLags',1,'Variance',garch(1,1),'Constant',0);
fit = estimate(Mdl,r) % fit the AR(1)-GARCH(1,1)
fit.Variance
[E0,V0] = infer(fit,r); % get the estimated errors and variances
%% Simulate periods ahead
[Y] = simulate(fit,2,'Y0',r,'E0',E0,'V0',V0, 'NumPaths',1e5);
ret = sum(Y); % compund the return
histfit(PV0*(exp(ret)-1),100,'normal')
title 'P&L distribution'
ret2d = prctile(ret,0.01); % get the 99% lowest return
VaR = PV0*(exp(ret2d)-1);
fprintf('Today portfolio value: %f\n2-days ahead 99%% VaR: %f\n',PV0,VaR);
Output:
ARIMA(1,0,0) Model:
--------------------
Conditional Probability Distribution: Gaussian
Standard t
Parameter Value Error Statistic
----------- ----------- ------------ -----------
Constant 0 Fixed Fixed
AR{1} -0.020272 0.0697731 -0.290541
GARCH(1,1) Conditional Variance Model:
----------------------------------------
Conditional Probability Distribution: Gaussian
Standard t
Parameter Value Error Statistic
----------- ----------- ------------ -----------
Constant 7.43521e-06 2.14546e-06 3.46556
GARCH{1} 0.653488 0.126818 5.15295
ARCH{1} 0.206016 0.0881823 2.33626
fit =
ARIMA(1,0,0) Model:
--------------------
Distribution: Name = 'Gaussian'
P: 1
D: 0
Q: 0
Constant: 0
AR: {-0.020272} at Lags [1]
SAR: {}
MA: {}
SMA: {}
Variance: [GARCH(1,1) Model]
ans =
GARCH(1,1) Conditional Variance Model:
--------------------------------------
Distribution: Name = 'Gaussian'
P: 1
Q: 1
Constant: 7.43521e-06
GARCH: {0.653488} at Lags [1]
ARCH: {0.206016} at Lags [1]
Today portfolio value: 2090.570000
2-days ahead 99% VaR: -79.986280
For a portfolio consisting of one share of S&P 500 index, I got the current value of \$2091, and 2-day VaR at 99% is \$80.
I also plotted the P&L distribution with Normal fit, so you can see that the normal distribution is not a very good fit. That's why you have to simulate when using GARCH.
You can also look at MATLAB's own GARCH example here.
The idea is that you fit AR(1)-GARCH(1,1) to the returns. Then you Monte-Carlo simulate the returns two days ahead. Then you compound two days of returns in each path, and find 0.01 percentile. This is how much or more you'll lose at 99% confidence in terms of returns, and in terms of dollars it's a simple arithmetic using the current portfolio value.
You can find worked out examples in this book: Carol Alexander, Market Risk Analysis, Volume IV, Value at Risk Models, February 2009
I know of at least five ways of initializing the volatility process:
1) Set it equal to $\varepsilon_{t-1}^2$,
2) The sample variance,
3) Unconditional variance of the model ($\alpha_0/(1-\alpha_1 - \alpha_2)$),
4) Allow it it to be an parameter to be estimated,
5) Backcasting with an exponential filter.
The topic is discussed in further detail here
Best Answer
To estimate $VaR$ of the one-step-ahead return (which is what you seem to need to do) you should do the following ($E_t$ denotes the conditional expectation given the values of the time series up to and including the moment $t$ have been observed):
Compute the forecast of the mean. In your case of $AR(p)$ process you should do the following: $$r_{t+1|t} = E_t(r_{t+1}) = E_t( \delta + \sum^p_{k=1} \phi_k r_{t+1-k} + \sigma_{t+1} \epsilon_{t+1})$$ which reduces to $r_{t+1|t} = \delta + \sum^p_{k=1} \phi_k r_{t+1-k}$, since $\epsilon_{t+1}$ has zero mean, and $\sigma^2_{t+1}$ is known given the observations up to $t$. In your case $$r_{t+1|t} = \delta + \phi_1 r_t$$
The volatility forecast is $\sigma^2_{t+1|t} = \alpha_0 + \alpha_1 \hat{z}^2_t + \beta \sigma^2_t$ where $\hat{z}_t = r_t - \delta - \phi_1 r_{t-1}$, and can be obtained similarly.
Now to get $VaR$ at the moment $t+1$, given the history up to $t$, you should just take the quantile of the standard normal, multiply it by the square root of $\sigma^2_{t+1|t}$ and then add $r_{t+1|t}$: $$VaR = r_{t+1|t} + \sqrt{\sigma^2_{t+1|t}} z_\alpha$$