\chapter{Численные методы}
\label{chap:numerical}

\section{Численная оптимизация}
\label{sec:hp-numopt}

В случае, когда эконометрист хочет использовать язык программирования,
такой как Hansl, а не просто полагаться на заранее подготовленные
подпрограммы, возникает необходимость применить какую-либо форму
численной оптимизации. Это может быть форма максимизации вероятности
или аналогичные методы индуктивной статистики. В другом варианте можно
использовать оптимизацию в более общем и абстрактном виде, например,
для решения задачи о выборе портфеля ценных бумаг или аналогичных
задач о распределении ресурсов.  Поскольку Hansl представляет собой
полный язык по Тьюрингу, в принципе любая числовая техника оптимизации
может программироваться в самом Hansl. Некоторые такие методы включены
в набор собственных функций программы, в интересах простоты ее
использования и эффективности. Они предназначены для решения наиболее
распространенных проблем, которые встречаются в экономике и
эконометрике, что неограниченно лишь оптимизацией дифференцируемых
функций.  В этой главе мы кратко рассмотрим, что конкретно предлагает
Hansl для решения задач следующего вида:
\[
\hat{\mathbf{x}} \equiv \argmax_{\mathbf{x} \in \Re^k} f(\mathbf{x}; \mathbf{a}),
\]
где $f(\mathbf{x}; \mathbf{a})$ --- функция от $\mathbf{x}$, форма
которой зависит от вектора параметров $\mathbf{a}$. Предполагается,
что целевая функция $f(\cdot)$ возвращает скаляру действительное
значение. В большинстве случаев также предполагается, что она
непрерывна и дифференцируема, хотя это необязательно. (Обратите
внимание, что хотя встроенные функции Hansl максимизируют заданную
целевую функцию, ее легко можно минимизировать, если просто
перевернуть знак $f(\cdot)$.)  Частный случай вышеизложенного
наступает, когда $\mathbf{x}$ --- это вектор параметров, а
$\mathbf{a}$ обозначает «данные». В этих случаях целевая функция
обычно представляет собой (логарифмическую) вероятность, а задача
заключается в ее оценке.  Для подобных случаев Hansl предлагает
несколько специальных конструкций, рассмотренных в
разделе~\ref{sec:est-blocks}. Здесь мы поговорим о более общих
проблемах. Тем не менее, важно отметить, что различия состоят только в
используемом синтаксисе Hansl; однако, математические алгоритмы,
которые использует gretl для решения задачи оптимизации, те же самые.
Читателю предлагается прочитать главу «Численные методы» Руководства
пользователя Gretl для полного представления о данной проблеме. Ниже
мы приведем лишь небольшой пример того, как это может выглядеть:

\begin{code}
function scalar Himmelblau(matrix x)
    /* extrema:
    f(3.0, 2.0) = 0.0, 
    f(-2.805118, 3.131312) = 0.0,
    f(-3.779310, -3.283186) = 0.0
    f(3.584428, -1.848126) = 0.0
    */
    scalar ret = (x[1]^2 + x[2] - 11)^2
    return -(ret + (x[1] + x[2]^2 - 7)^2)
end function

# ----------------------------------------------------------------------

set max_verbose 1

matrix theta1 = { 0, 0 }
y1 = BFGSmax(theta1, "Himmelblau(theta1)")
matrix theta2 = { 0, -1 }
y2 = NRmax(theta2, "Himmelblau(theta2)")

print y1 y2 theta1 theta2
\end{code}
Мы используем для иллюстрации классическую «неприятную» функцию из
литературы по численной оптимизации, а именно функцию Химмельблау,
которая имеет четыре различных минимума
$f(x, y) = (x^2+y-11)^2 + (x+y^2-7)^2$. Алгоритм выглядит следующим
образом:

\begin{enumerate}
\item Сначала мы определяем функцию для оптимизации: она должна
  возвращать скалярное значение и иметь среди своих аргументов вектор
  для оптимизации. В данном конкретном случае это единственный
  аргумент, но могли быть также и другие. Поскольку в этом случае мы
  минимизируем функцию, программа возвращает отрицательное значение
  функции Himmelblau.
\item Затем мы задаем \verb|max_verbose| как 1. Это еще один пример
  использования команды \cmd{set} ; сама команда предлагает: «давайте
  посмотрим, как проходят итерации», и по умолчанию она равна 0. С
  помощью команды \cmd{set} с соответствующими параметрами, можно
  управлять несколькими функциями процесса оптимизации, такими как
  числовые значения допусков, визуализация итераций и т. д.
\item Определяем $\theta_1 = [0, 0]$ как начальную точку.
\item Применяем функцию \cmd{BFGSmax}; она будет искать максимум с
  помощью техники BFGS. Ее базовый синтаксис выглядит так:
  \texttt{BFGSmax(arg1, arg2)}, где \texttt{arg1} --- вектор,
  содержащий переменную оптимизации, а \texttt{arg2} --- это строка,
  содержащая вызов функции, которую нужно максимизировать. BFGS
  подставляет значения $\theta_1$, пока максимум не будет
  достигнут. При успешном завершении вектор \texttt{theta1} будет
  содержать конечную точку. (Примечание: об этом сказано гораздо
  больше в Руководстве пользователя Gretl и Справочнике по командам
  Gretl.)
\item Затем мы решаем ту же проблему, но с другой отправной точкой и
  другой техникой оптимизации. Мы начинаем с $\theta_2 = [0, -1]$ и
  используем метод Ньютона – Рафсона вместо BFGS, применяя
  \cmd{NRmax()} вместо \cmd{BFGSmax()}. Синтаксис тот же.
\item Выводим результаты.
\end{enumerate}

Данный вывод содержится в ~\ref{tab:optim-output}. Обратите внимание,
что второй раз получается иной локальный оптимум, чем в первый
раз. Это следствие того, что мы начали тот же алгоритм с другой
отправной точки. В этом примере использовались числовые производные,
но вы можете использовать аналитически вычисленные производные для
обоих методов, если у вас есть функция Hansl для них; см. Руководство
пользователя Gretl для получения более подробной информации.

\begin{table}[ht]
  \begin{footnotesize}
\begin{scode}
? matrix theta1 = { 0, 0 }
Replaced matrix theta1
? y1 = BFGSmax(theta1, "Himmelblau(11, theta1)")
Iteration 1: Criterion = -170.000000000
Parameters:       0.0000      0.0000
Gradients:        14.000      22.000 (norm 0.00e+00)

Iteration 2: Criterion = -128.264504038 (steplength = 0.04)
Parameters:      0.56000     0.88000
Gradients:        33.298      39.556 (norm 5.17e+00)

...

--- FINAL VALUES: 
Criterion = -1.83015730011e-28 (steplength = 0.0016)
Parameters:       3.0000      2.0000
Gradients:    1.7231e-13 -3.7481e-13 (norm 7.96e-07)

Function evaluations: 39
Evaluations of gradient: 16
Replaced scalar y1 = -1.83016e-28
? matrix theta2 = { 0, -1 }
Replaced matrix theta2
? y2 = NRmax(theta2, "Himmelblau(11, theta2)")
Iteration 1: Criterion = -179.999876556 (steplength = 1)
Parameters:   1.0287e-05     -1.0000
Gradients:        12.000  2.8422e-06 (norm 7.95e-03)

Iteration 2: Criterion = -175.440691085 (steplength = 1)
Parameters:      0.25534     -1.0000
Gradients:        12.000  4.5475e-05 (norm 1.24e+00)

...

--- FINAL VALUES: 
Criterion = -3.77420797114e-22 (steplength = 1)
Parameters:       3.5844     -1.8481
Gradients:   -2.6649e-10  2.9536e-11 (norm 2.25e-05)

Gradient within tolerance (1e-07)
Replaced scalar y2 = -1.05814e-07
? print y1 y2 theta1 theta2

             y1 = -1.8301573e-28

             y2 = -1.0581385e-07

theta1 (1 x 2)

  3   2 

theta2 (1 x 2)

      3.5844      -1.8481 
\end{scode}
    
  \end{footnotesize}
  \caption{Результат максимизации}
  \label{tab:optim-output}
\end{table}
Язык Hansl предлагает вашему вниманию следующие методы оптимизации:
\begin{itemize}
\item BFGS через функцию \cmd{BFGSmax()}. В большинстве случаев это
  лучший компромисс между производительностью и
  надежностью. Предполагается, что функция максимизации
  дифференцируема и будет аппроксимировать кривую за счет
  использования изменений градиента между итерациями. Вы можете
  снабдить ее аналитически вычисленным градиентом для большей скорости
  и точности, но если его нет, то первые производные будут вычисляться
  численно.
\item Метод Ньютона --- Рафсона с помощью функции \cmd{NRmax()}. На
  самом деле это название немного вводит в заблуждение. Метод должен
  был называться чем-то вроде «на основе кривизны функции», поскольку
  он основан на итерациях
  \[
    x_{i+1} = -\lambda_i C(x_i)^{-1} g(x_i)
  \]
  где $g(x)$ градиент, а $C(x_i)$ некоторая мера кривизны функции для
  оптимизации; если $C(x)$ --- матрица Гессиана, то вы получите
  функцию Ньютона – Рафсона. Опять же, вы можете закодировать свои
  собственные функции для $g(\cdot)$ и $C(\cdot)$, но если вы этого не
  сделаете, то будут использоваться соответственно численные
  приближения градиента и гессиана. Другие популярные методы
  оптимизации, такие как BHHH и алгоритм подсчета очков, могут быть
  применены путем добавления в \cmd{NRmax()} соответствующей матрицы
  кривизны $C(\cdot)$. Этот метод очень эффективен, когда он работает,
  но довольно хрупок в использовании: например, если $C(x_i)$
  оказывается неотрицательно определенной на некоторой итерации,
  сходимость может стать проблематичной.
\item Методы без производных: единственный метод, который в настоящее
  время предлагает Hansl, --- это алгоритм имитации отжига (simulated
  annealing) через функцию \cmd{simann()}, но реализация алгоритма
  Нелдера – Мида (также известного как метод «амебы») лишь вопрос
  времени. Эти методы работают даже тогда, когда функция максимизации
  имеет некоторую форму разрыва или не везде дифференцируема; однако
  они могут работать очень медленно и перегружать процессор.
\end{itemize}

\section{Численное дифференцирование}
\label{sec:hp-numdiff}

Для численного дифференцирования возможно применить функцию
\texttt{fdjac}. Например:

\begin{code}
set echo off
set messages off

function scalar beta(scalar x, scalar a, scalar b)
    return x^(a-1) * (1-x)^(b-1)
end function

function scalar ad_beta(scalar x, scalar a, scalar b)
    scalar g = beta(x, a-1, b-1)
    f1 = (a-1) * (1-x)
    f2 = (b-1) * x
    return (f1 - f2) * g
end function

function scalar nd_beta(scalar x, scalar a, scalar b)
    matrix mx = {x}
    return fdjac(mx, beta(mx, a, b))
end function

a = 3.5
b = 2.5

loop for (x=0; x<=1; x+=0.1)
    printf "x = %3.1f; beta(x) = %7.5f, ", x, beta(x, a, b)
    A = ad_beta(x, a, b)
    N = nd_beta(x, a, b)
    printf "analytical der. = %8.5f, numerical der. = %8.5f\n", A, N
endloop
\end{code}

возвращает
\begin{code}
x = 0.0; beta(x) = 0.00000, analytical der. =  0.00000, numerical der. =  0.00000
x = 0.1; beta(x) = 0.00270, analytical der. =  0.06300, numerical der. =  0.06300
x = 0.2; beta(x) = 0.01280, analytical der. =  0.13600, numerical der. =  0.13600
x = 0.3; beta(x) = 0.02887, analytical der. =  0.17872, numerical der. =  0.17872
x = 0.4; beta(x) = 0.04703, analytical der. =  0.17636, numerical der. =  0.17636
x = 0.5; beta(x) = 0.06250, analytical der. =  0.12500, numerical der. =  0.12500
x = 0.6; beta(x) = 0.07055, analytical der. =  0.02939, numerical der. =  0.02939
x = 0.7; beta(x) = 0.06736, analytical der. = -0.09623, numerical der. = -0.09623
x = 0.8; beta(x) = 0.05120, analytical der. = -0.22400, numerical der. = -0.22400
x = 0.9; beta(x) = 0.02430, analytical der. = -0.29700, numerical der. = -0.29700
x = 1.0; beta(x) = 0.00000, analytical der. = -0.00000, numerical der. =       NA
\end{code}

Подробности об используемом алгоритме можно найти в Справочнике по
командам Gretl. Достаточно сказать, что здесь параметр
\texttt{fdjac\_quality} задан от 0 до 2. Значение по умолчанию --- 0,
который дает приближение прямой разницы: это самый быстрый алгоритм,
но иногда он может быть недостаточно точным. Значение 1 дает
двустороннюю разницу, а 2 использует экстраполяцию Ричардсона. По мере
того, как значение параметра повышается, растет и точность, но метод
значительно больше загружает процессор.

% \section{Random number generation}

% \begin{itemize}
% \item Mersenne Twister in its various incarnations
% \item Ziggurat vs Box--Muller
% \item Other distributions
% \end{itemize}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "hansl-primer"
%%% End: 

