SAS: работа с матрицами с использованием IML

Итак, я пытаюсь рассчитать эту формулу, но результаты странные. Элементы очень большие, поэтому я не уверен, где я ошибся. Я приложил фото формулы: введите здесь описание изображения

и вот мой код:

*calculating mu_sum and sigma_sum;
T_hat=180;
mu_sum_first_part={0,0,0,0};
mu_sum_second_part={0,0,0,0};
mu_sum={0,0,0,0};

*calculating mu_sum;
do i = 0 to T_hat;
    term=(T_hat - i)*(B0**i)*a;
    mu_sum_first_part = mu_sum_first_part + term;
end;
do i=1 to T_hat; 
    term =B0**i;
    mu_sum_second_part = mu_sum_second_part + term;
end;
mu_sum = mu_sum_first_part + mu_sum_second_part*zt;
print mu_sum;

*calculating sigma_sum;
term=I(4);
sigma_sum=sigma;
do j=1 to T_hat;
    term = term + B0**j;
    sigma_sum = sigma_sum + (term*sigma*(term`));
end;
print sigma_sum;

Я знаю, что это долго, но, пожалуйста, помогите!


person duckman    schedule 12.07.2017    source источник
comment
Что это за расчет, хотелось бы посмотреть? Является ли Bo скаляром, вектором или матрицей? Я бы подумал о векторе, но I+Bo в этом смысле не имеет смысла.   -  person DomPazz    schedule 12.07.2017
comment
@DomPazz Я выполняю задачу по распределению активов. это из статьи Барбериса (JF, 2000) Инвестирование в долгосрочной перспективе, когда доходность предсказуема. Bo представляет собой матрицу 4x4, в которой первый столбец содержит нули. альфа — это матрица 4x1, а сигма — матрица 4x4.   -  person duckman    schedule 13.07.2017
comment
Я не тянул бумагу, является ли оператор степени на матрице степенью матрицы или поэлементной мощностью? IE B^2 = B`*B или B#B, где # — поэлементное умножение   -  person DomPazz    schedule 13.07.2017
comment
Оператор двойной звезды представляет собой многократное умножение матриц.   -  person Rick    schedule 13.07.2017
comment
@DomPazz оператор мощности - это матричная мощность, а не поэлементная мощность   -  person duckman    schedule 15.07.2017


Ответы (2)


Первое, что бросается в глаза, это то, что первый член вашего цикла в mu имеет 1 слишком много:

do i = 0 to T_hat;
    term=(T_hat - i)*(B0**i)*a;
    mu_sum_first_part = mu_sum_first_part + term;
end;

Должно быть:

do i = 0 to T_hat-1;
    term=(T_hat - i)*(B0**i)*a;
    mu_sum_first_part = mu_sum_first_part + term;
end;
person DomPazz    schedule 13.07.2017
comment
оба должны быть такими же, как (T-hat - i) = 0, когда i = T_hat, поэтому последний член равен 0. - person duckman; 15.07.2017

В вашей программе нет ничего математически неправильного. Когда вы возводите матрицу в 180-ю степень, вы не должны удивляться, увидев очень большие или очень маленькие значения. Например, если вы позволите B0 = {0 1 0 0, 0 0 1 0, 0 0 0 1, 0 1 1 1 }; тогда элементами B0**T являются O( 1E47 ). Если вы разделите B0 на 2 и возведете результат в 180-ю степень, то элементы будут O(1E-8).

Предположительно, эти формулы предназначены для матриц B0, имеющих специальную структуру, например ||B0**n|| --> 0 как n --> бесконечность. В противном случае степенной ряд не сойдется. Я предлагаю вам перепроверить, что B0, который вы используете, удовлетворяет предположениям ссылки.

Вы не спрашивали об эффективности, но было бы разумно вычислить усеченный ряд мощности, используя метод Хорнера в SAS/IML вместо явного формирования степеней B0.

person Rick    schedule 13.07.2017