Matlab: AR(1) и GARCH(1,1) своими руками.

by

Часть 0. AR(1)

Генератор AR(1)-процесса. Оценка параметров методом максимального правдоподобия.


function [y] = generateAR1(T, c, fi, sigma2)
 % Generates AR(1) – time series process
 % T = number of observations
 % с = constant of the process
 % fi = coefficient
 % sigma2 = variance of epsilon; epsilon is a noise of the prosess
 epsilon = (sigma2^0.5).*randn(T,1);
 y = zeros(T,1);
 y(1) = ((sigma2/(1-fi^2))^0.5).*randn(1,1) + c/(1-fi);
 for t = 2:T
 y(t) = c + fi*y(t-1) + epsilon(t);
 end

y(t) = 0.8 * y(t-1) + e(t)


function [theta_ml, LogLikelihoodF] = LLestimateAR1 (y)
 % returns ML - estimators and value of log-likelihood function

fun1 = @(theta) LLAR1(y, theta);
 theta0 = [1 0.5 1]';
 [theta_ml, negative_LogLikelihoodF] = fminsearch(fun1, theta0);
 LogLikelihoodF = -negative_LogLikelihoodF;

function [negative_LL] = LLAR1(y, theta)
 T = max(size(y));
 c = theta(1);
 fi = theta(2);
 sigma2=theta(3);

LL1 = 0.5*log(2*pi);
 LL2 = 0.5*log(sigma2/(1-fi^2));
 LL3 = ((y(1)-(c/(1-fi)))^2)*(2*sigma2/(1-fi^2))^(-1);
 LL4 = ((T-1)/2)*log(2*pi);
 LL5 = ((T-1)/2)*log(sigma2);
 LL6 = 0;
 for t=2:T
 LL6 = LL6 + ((y(t)-c-fi*y(t-1))^2)/(2*sigma2);
 end
negative_LL = LL1+LL2+LL3+LL4+LL5+LL6;

Если сначала сгенерируете AR(1)-процесс.

>> [y] = generateAR1(500, 0, 0.8, 1);

А потом оцените параметры.

>> [theta_ml, LogLikelihoodF] = LLestimateAR1 (y);

То получите нечто вроде этого:

theta_ml =

0.0209
0.7990
0.9890

LogLikelihoodF = -707.2230

Часть 1. Генератор GARCH(1,1)-процесса.


function [y, sigma] = generate_GARCH11(N, c, w, delta, gamma)
y     = zeros(N, 1);
epsilon   = zeros(N, 1);
sigma = zeros(N, 1);
epsilon(1) = sqrt(w/(1 - delta - gamma)) * randn;
sigma(1) = sqrt(w/(1 - delta - gamma));
ksi = randn(N, 1);

for i=2:N
sigma(i) = (w + delta*sigma(i-1)^2 + gamma*epsilon(i-1)^2)^0.5;
epsilon(i) = sigma(i)*ksi(i);
y(i) = c + epsilon(i);
end

Функция, которая оценивает параметры GARCH(1,1)-процесса.


function [c_ml, w_ml, delta_ml, gamma_ml, LLvalue] = estimate_GARCH11_fmincon(y)

fun1 = @(theta) garchml(theta, y);

%__________________________________________________________________________
% Начальная точка в процессе оптимизации.
c0 = 1;
w0 = 1;
delta0 = 0.2;
gamma0 = 0.2;
theta0 = [c0 w0 delta0 gamma0];
%__________________________________________________________________________

%__________________________________________________________________________
H = 10^4; % Верхнее ограничение.
h = 10^(-4); % Нижнее ограничение.
%__________________________________________________________________________

%__________________________________________________________________________
% Параметры для функции fmincon.
fun = fun1;
x0 = theta0;
A = [0 0 1 1];
b = 1 - h;
Aeq = [];
beq = [];
lb = [-H h h h]';
ub = [H H H H]';
nonlcon = [];
options = optimset('Algorithm', 'interior-point', 'Display', 'off');
%__________________________________________________________________________

%__________________________________________________________________________
% Минимизация логарифмической функции правдоподобия при помощи функции
% fmincon.
[theta_ml, negLLvalue] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
%__________________________________________________________________________

%__________________________________________________________________________
c_ml = theta_ml(1);
w_ml = theta_ml(2);
delta_ml = theta_ml(3);
gamma_ml = theta_ml(4);
LLvalue = -negLLvalue;
%__________________________________________________________________________

function f = garchml(teta, y)

c = teta(1);
w = teta(2);
delta = teta(3);
gamma = teta(4);

N = max(size(y));
ym = sum(y)/N;

eps = zeros(N,1);
sigma = zeros(N,1);
sigma(1) = (sum((y-ym).^2)/N)^0.5;

for i=2:N
eps(i) = y(i) - c;
sigma(i) = (w + delta*sigma(i-1)^2 + gamma*eps(i-1)^2)^0.5;
end

f = sum(log(sigma.^2)) + sum((eps.^2)./(sigma.^2));

Результаты следующие.

Расчет № 1.

>> [y, sigma] = generate_GARCH11(1000, 0, 1, 0.5, 0.3);

c_ml = -0.0291

w_ml = 1.0066

delta_ml = 0.4782

gamma_ml = 0.3239

LLvalue = -2.4270e+003

Расчет № 2.
>> [y, sigma] = generate_GARCH11(100000, 0, 1, 0.5, 0.3);
>> [c_ml, w_ml, delta_ml, gamma_ml, LLvalue] = estimate_GARCH11_fmincon(y)

c_ml = 0.0015

w_ml = 0.9729

delta_ml = 0.5084

gamma_ml = 0.2966

LLvalue = -2.4383e+005

Часть 2.
Функция, при помощи которой оцениваются коэффициенты и ковариационная матрица коэффициентов.

Данная функция написана с использованием материала, который любезно предоставил в свободный доступ А.А. Цыплаков.

http://www.nsu.ru/ef/tsy/ecmr/garch/short/GARCH.pdf


function [theta_ml_tsy, theta_ml_tsy_std, cov_matr_theta, iter] = estimate_GARCH11_tsy_ver3(y)
iter_max = 100;
iter = 1;

T = length(y);
c_0 = mean(y);
uncond_sigma2 = sum((y - mean(y)).^2)/T;
delta_0 = 1/4;
gamma_0 = 1/4;
omega_0 = uncond_sigma2 * (1 - delta_0 - gamma_0);
theta_0 = [c_0, omega_0, delta_0, gamma_0]';
k = length(theta_0);

epsilon = y - c_0;
sigma2 = zeros(T, 1);
sigma2(1,1) = omega_0 + delta_0*uncond_sigma2 + gamma_0*uncond_sigma2;
for t = 2:T
sigma2(t,1) = omega_0 + delta_0*sigma2(t-1,1) + gamma_0*(epsilon(t-1,1)^2);
end

grad_epsilon = zeros(T,k);
for t = 1:T
grad_epsilon(t,:) = [-1, 0, 0, 0];
end

grad_sigma2 = zeros(T,k);
grad_sigma2(1,:) = [0, 1, uncond_sigma2, uncond_sigma2];
for t = 2:T
grad_sigma2(t,:) = [0, 1, sigma2(t-1,1), epsilon(t-1,1)^2] + delta_0*grad_sigma2(t-1,:) + gamma_0*epsilon(t-1,1)*grad_epsilon(t-1,:);
end

xi = epsilon./sqrt(sigma2);
nu = xi.^2 - 1;
Q = -grad_epsilon./[sqrt(sigma2), sqrt(sigma2), sqrt(sigma2), sqrt(sigma2)];
S = grad_sigma2./[sigma2, sigma2, sigma2, sigma2];

delta_theta = (Q'*Q + 0.5*(S')*S)\(Q'*xi + 0.5*S'*nu);

Y_secondary = [xi; (1/sqrt(2))*nu];
X_secondary = [Q; (1/sqrt(2))*S];
TSS_secondary = sum((Y_secondary - mean(Y_secondary)).^2);
RSS_secondary = (Y_secondary - X_secondary*delta_theta)'*(Y_secondary - X_secondary*delta_theta);
Rsq_secondary = 1 - (RSS_secondary/TSS_secondary);

stop = (Rsq_secondary < 10^(-5));

if (stop == 1)
theta_ml_tsy = theta_0;
cov_matr_theta = (RSS_secondary/(2*T))*inv(Q'*Q + 0.5*(S')*S);
% cov_matr_theta = (TSS_secondary/(2*T))*inv(Q'*Q + 0.5*(S')*S);
theta_ml_tsy_std = [sqrt(cov_matr_theta(1,1)), sqrt(cov_matr_theta(2,2)), sqrt(cov_matr_theta(3,3)), sqrt(cov_matr_theta(4,4))]';
return
end

lambda = 1;
theta_temp = theta_0 + lambda*delta_theta;
ok = (GARCH11_LF(theta_temp, y) > GARCH11_LF(theta_0, y))&&(isinG(theta_temp));
while (ok == 0)
lambda = 0.9*lambda;
theta_temp = theta_0 + lambda*delta_theta;
ok = (GARCH11_LF(theta_temp, y) > GARCH11_LF(theta_0, y))&&(isinG(theta_temp));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:iter_max
iter = iter + 1;
c_0 = theta_temp(1);
omega_0 = theta_temp(2);
delta_0 = theta_temp(3);
gamma_0 = theta_temp(4);
theta_0 = [c_0, omega_0, delta_0, gamma_0]';
k = length(theta_0);

epsilon = y - c_0;
sigma2 = zeros(T, 1);
sigma2(1,1) = omega_0 + delta_0*uncond_sigma2 + gamma_0*uncond_sigma2;
for t = 2:T
sigma2(t,1) = omega_0 + delta_0*sigma2(t-1,1) + gamma_0*(epsilon(t-1,1)^2);
end

grad_epsilon = zeros(T,k);
for t = 1:T
grad_epsilon(t,:) = [-1, 0, 0, 0];
end

grad_sigma2 = zeros(T,k);
grad_sigma2(1,:) = [0, 1, uncond_sigma2, uncond_sigma2];
for t = 2:T
grad_sigma2(t,:) = [0, 1, sigma2(t-1,1), epsilon(t-1,1)^2] + delta_0*grad_sigma2(t-1,:) + gamma_0*epsilon(t-1,1)*grad_epsilon(t-1,:);
end

xi = epsilon./sqrt(sigma2);
nu = xi.^2 - 1;
Q = -grad_epsilon./[sqrt(sigma2), sqrt(sigma2), sqrt(sigma2), sqrt(sigma2)];
S = grad_sigma2./[sigma2, sigma2, sigma2, sigma2];

delta_theta = (Q'*Q + 0.5*(S')*S)\(Q'*xi + 0.5*S'*nu);

Y_secondary = [xi; (1/sqrt(2))*nu];
X_secondary = [Q; (1/sqrt(2))*S];
TSS_secondary = sum((Y_secondary - mean(Y_secondary)).^2);
RSS_secondary = (Y_secondary - X_secondary*delta_theta)'*(Y_secondary - X_secondary*delta_theta);
Rsq_secondary = 1 - (RSS_secondary/TSS_secondary);

stop = (Rsq_secondary < 10^(-5));

if (stop == 1)
theta_ml_tsy = theta_0;
cov_matr_theta = (RSS_secondary/(2*T))*inv(Q'*Q + 0.5*(S')*S);
% cov_matr_theta = (TSS_secondary/(2*T))*inv(Q'*Q + 0.5*(S')*S);
theta_ml_tsy_std = [sqrt(cov_matr_theta(1,1)), sqrt(cov_matr_theta(2,2)), sqrt(cov_matr_theta(3,3)), sqrt(cov_matr_theta(4,4))]';
return
end

lambda = 1;
theta_temp = theta_0 + lambda*delta_theta;
ok = (GARCH11_LF(theta_temp, y) > GARCH11_LF(theta_0, y))&&(isinG(theta_temp));
while (ok == 0)
lambda = 0.9*lambda;
theta_temp = theta_0 + lambda*delta_theta;
ok = (GARCH11_LF(theta_temp, y) > GARCH11_LF(theta_0, y))&&(isinG(theta_temp));
end

end

Замечание.

В текущую папку, где находится эта программа поместите следующую функцию.


function [rez] = isinG(x)
h = 10^(-5);
H = 10^2;
if (x(2) > h)&&(x(2) < H)&&(x(3) > h)&&(x(4) > h)&&(x(3) + x(4) < 1-h)
rez = 1;
else rez = 0;
end

Пример расчета.

>> [y, sigma] = generate_GARCH11(1000, 0, 1, 0.5, 0.3);
>> [theta_ml_tsy, cov_matr_theta, iter] = estimate_GARCH11_tsy_ver2(y)

theta_ml_tsy =

0.0562
1.2541
0.4114
0.3783

cov_matr_theta =

0.0077    0.0000    0.0000   -0.0000
0.0000    0.1089   -0.0258    0.0071
0.0000   -0.0258    0.0088   -0.0047
-0.0000    0.0071   -0.0047    0.0059

iter = 4

Реклама

Метки: ,

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход / Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход / Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход / Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход / Изменить )

Connecting to %s


%d такие блоггеры, как: