Un exercici d’optimització amb solvers no lineals de LibreOffice Calc i Octave

Qui segueixi aquest bloc ja sabrà que soc un “aficionat” als MOOC de Coursera. Doncs bé, recentment n’he acabat un de “Financial Engineering and Risk Management I” (l”I’ és perquè tindrà continuació aquest gener).

El curs ha tingut un important contingut pràctic que ha demanat l’utilització del full de càlcul.

Excel, com no, ha estat el full de càlcul recomanat i els fulls d’exercicis i plantilles del material didàctic han estat servits en format .xlsx. També es suggeria des de l’staff del curs la realització dels exercicis amb matlab, o amb Python, o amb les alternatives de programari lliure: OpenOffice.org/LibreOffice Calc, o gnumeric, o Octave.

Personalment he seguit el curs amb el LibreOffice Calc. No hi ha hagut cap problema i el full de càlcul ha permès la resolució de tots els problemes. Ara bé, sí que he trobat diferències de rendiment en la resolució d’un problema de programació no-lineal que m’ha cridat l’atenció.

El problema

El problema que es plantejava tenia una primera part que era la construcció i calibració d’un model de Black-Derman-Toy de interès ri,j de 10 períodes.

El model pretén fer una estimació de les possibles evolucions del tipus d’interès en uns períodes de temps. Després aquest tipus d’interès es fa servir per estimar l’evoloció dels preus de bons del tresor per a, finalment, calcular el preu de productes financers derivats.

L’esquema anterior es tradueix en tres blocs al full de càlcul, un primer bloc que correspon a l’estimació de l’evolució del tipus d’interès,

Un segon bloc que, fent servir l’estimació anterior, fa l’estimació dels preus dels bons.

I un tercer bloc que partint del bloc del preu dels bons, calcula el preu del derivat.

Per a que el sistema funcioni cal ajustar els dos primers blocs per a que el valor estimat dels preus dels bons sigui coherent amb els tipus d’interès observats als mercat.

En el model BDT els tipus d’interès segueixen una llei exponencial r = a * exp(b * t). Es fixa “b” per a tot el bloc, que segueix una estructura de diagrama de Lattice; i “a” és variable i depèn del període (cada columna correspon a un període).

la següent imatge, amb els rastres en blau de l’eina detectiu del Calc que permet marcar les precedències de les fórmules, il·lustren l’explicació que ve a continuació:

escriptori 1_004

Es parteix d’una estimació inicial de valors “a”. En el cas de l’exercici, és la fila en groc. Els valors d'”a” observats al mercat són a la fila en verd.

El bloc de tipus d’interès es construeix simplement aplicant la fórmula r = a * exp(b * t) pera cada període, però donant-li aquesta estructura triangular que serà útil per a la construcció del següent bloc.

El bloc de preu dels bons es construeix partint d’un preu unitari en temps = 0 i avançant cap a t+1 mitjançant les forward equations, que, al seu cop, fan servir els tipus d’interès del període t. Aquesta és l’explicació de l’estructura de tipus d’interès que he construït abans. Al final tenim uns estats que depenen del període t, i de la història i les variacions prèvies dels tipus d’interès.

Tots aquests valors probables del bon per a un període permeten per addició obtenir un valor estimat del preu del bon per a aquell període. Hem calculat, doncs, un preu del bon basat en l’estimació dels tipus d’interès. És la fila en marró.

Però, de fet, el que ens interessa és que el model s’ajusti als valors observats al mercat. Per tant, el que es fa és calcular el tipus d’interès estimat per al període que es dedueix del preu estimat del bon. És la primera fila en blau.

El que volem és minimitzar la diferència entre el valor estimat i el de mercat. Per això és calcula el quadrat de la diferència de cada període (darrera línea en blau) i se sumen tots els períodes per a obtenir una cel·la objectiu, en vermell, que vull minimitzar.

Els solvers

El problema plantejat és no lineal. Aquesta mena de problemes es poden resoldre amb els solvers no lineals que venen amb alguns fulls de càlcul. No amb tots. Per exemple, el Google Calc no té un solver no lineal (Sí que el té lineal). I les versions antigues del OpenOffice.org tampoc el tenien “de fàbrica” i cali afegir un plugin de Sun.

Però Excel, Gnumeric i Calc de LibreOffice/OpenOffice.org sí que munten solvers no lineals, tot i que en el cas d’Excel i de Calc, cal activar-los explícitament.

El funcionament dels tres solvers és similar: es tracta de dir quina és la cel·la objectiu; indicar quina mena d’optimització es vol fer: maximitzar, minimitzar o igualar a un valor (gnumeric no permet aquesta última opció); indicar quines cel·les cal variar per obtenir el resultat; i proporcionar un valor inicial per arrencar les iteracions. A més tots tres solvers permeten afegir condicions addicionals (encara que sembla que el Calc no acaba de fer-ne cas d’aquestes condicions). En tots tres casos es possible definir el grau de precisió de la resposta; que fer si l’algorisme s’encalla o posar un límit en el nombre d’iteracions.

La resolució de problemes de programació no lineal es pot atacar amb diferents algorismes. La diferència entre els diferents solvers és, precisament, l’algorisme que fan servir

L’algorisme de l’Excel resol el problema amb gran rapidesa. El de Gnumeric és el que segueix en velocitat de resolució. Libre Office Calc permet triar dos algorismes diferents i permet fer un seguiment detallat del procés. Si simplement es vol obtenir el resultat i prenent els valors per defecte del solver no he notat una diferència sensible de temps entre els dos algorismes proposats, però sí que he notat que trigava més a obtenir resultats comparat amb els altres fulls de càlcul.

Caldria veure si el problema és d’eficiència de l’algorisme, o deficiència (sense apòstrof) de la implementació.

En tot cas, fent servir l'”SCO Evolutionary Algorithm” i ajustant algun paràmetre de precisió, obtinc el resultat:

escriptori 1_006

La solució amb octave

En comptes d’encarar el problema amb el full de càlcul també es pot fer solucionar amb llenguatges de programació. Segurament el llenguatge estrella per a aquesta mena de problemes deu ser Matlab. O la seva alternativa de programari lliure, Octave.

Com alternativa, també es pot fer servir un llenguatge general més una llibreria. L’aproximació més comú avui segurament seria la combinació de Python amb llibreria SciPy/NumPy; o Java, també amb alguna llibreria que proporcioni el solver. Per descomptat, sempre es pot intentar l’opció DIY -do it yourself- que segur que també és molt interessant. A Gnumeric estan buscant algú que s’animi a millorar el seu solver 😉

La meva solució al problema es basa en la funció sqp. En la seva forma d’utilització més senzilla, se li passa una funció a minimitzar, en el meu cas bdt_objective i un vector amb una aproximació a la solució.

sqp admet més paràmetres. Aquesta és la pàgina de documentació de la funció que reprodueixo a continuació

25.3 Nonlinear Programming

Octave can also perform general nonlinear minimization using a successive quadratic programming solver.


— Function File: [x, obj, info, iter, nf, lambda] = sqp (x0, phi)
— Function File: [...] = sqp (x0, phi, g)
— Function File: [...] = sqp (x0, phi, g, h)
— Function File: [...] = sqp (x0, phi, g, h, lb, ub)
— Function File: [...] = sqp (x0, phi, g, h, lb, ub, maxiter)
— Function File: [...] = sqp (x0, phi, g, h, lb, ub, maxiter, tol)
Solve the nonlinear program

          min phi (x)
           x
subject to

          g(x)  = 0
          h(x) >= 0
          lb <= x <= ub
using a sequential quadratic programming method.

The first argument is the initial guess for the vector x0.

The second argument is a function handle pointing to the objective function phi. The objective function must accept one vector argument and return a scalar.

The second argument may also be a 2- or 3-element cell array of function handles. The first element should point to the objective function, the second should point to a function that computes the gradient of the objective function, and the third should point to a function that computes the Hessian of the objective function. If the gradient function is not supplied, the gradient is computed by finite differences. If the Hessian function is not supplied, a BFGS update formula is used to approximate the Hessian.

When supplied, the gradient function phi{2} must accept one vector argument and return a vector. When supplied, the Hessian function phi{3} must accept one vector argument and return a matrix.

The third and fourth arguments g and h are function handles pointing to functions that compute the equality constraints and the inequality constraints, respectively. If the problem does not have equality (or inequality) constraints, then use an empty matrix ([]) for g (or h). When supplied, these equality and inequality constraint functions must accept one vector argument and return a vector.

The third and fourth arguments may also be 2-element cell arrays of function handles. The first element should point to the constraint function and the second should point to a function that computes the gradient of the constraint function:

                      [ d f(x)   d f(x)        d f(x) ]
          transpose ( [ ------   -----   ...   ------ ] )
                      [  dx_1     dx_2          dx_N  ]
The fifth and sixth arguments, lb and ub, contain lower and upper bounds on x. These must be consistent with the equality and inequality constraints g and h. If the arguments are vectors then x(i) is bound by lb(i) and ub(i). A bound can also be a scalar in which case all elements of x will share the same bound. If only one bound (lb, ub) is specified then the other will default to (-realmax, +realmax).

The seventh argument maxiter specifies the maximum number of iterations. The default value is 100.

The eighth argument tol specifies the tolerance for the stopping criteria. The default value is sqrt(eps).

The value returned in info may be one of the following:

101
The algorithm terminated normally. Either all constraints meet the requested tolerance, or the stepsize, delta x, is less than tol * norm (x). 
102
The BFGS update failed. 
103
The maximum number of iterations was reached.
An example of calling sqp:

          function r = g (x)
            r = [ sumsq(x)-10;
                  x(2)*x(3)-5*x(4)*x(5);
                  x(1)^3+x(2)^3+1 ];
          endfunction
          
          function obj = phi (x)
            obj = exp (prod (x)) - 0.5*(x(1)^3+x(2)^3+1)^2;
          endfunction
          
          x0 = [-1.8; 1.7; 1.9; -0.8; -0.8];
          
          [x, obj, info, iter, nf, lambda] = sqp (x0, @phi, @g, [])
          
          x =
          
            -1.71714
             1.59571
             1.82725
            -0.76364
            -0.76364
          
          obj = 0.053950
          info = 101
          iter = 8
          nf = 10
          lambda =
          
            -0.0401627
             0.0379578
            -0.0052227
See also: qp.

La dificultat, doncs, es redueix a implementar correctament la funció objectiu. Crec que el codi es prou entenedor. El procés és similar al seguit amb el Calc i les imatges prèvies poden ajudar a explicar el codi.

function objective = bdt_objective(a)
  spot_rate = [3, 3.1, 3.2, 3.3, 3.4, 3.5, 3.55, 3.6, 3.65, 3.7];
  b = 0.05;
  q = 0.5;

  short_rate_lattice = zeros(10,10);
  bond_price_lattice = zeros(11,11);
  bdt_model_zcb_prices = [];
  bdt_model_spot_rates = [];
  squared_diffs = [];
  objective = 0;

  for row = [1:10]
    for col = [1:10]
      if (col + row >= 11)
        short_rate_lattice(row, col) = a(col)*exp(b*(10 - row));
      else
        short_rate_lattice(row, col) = 0;
      endif
    endfor 
  endfor

  bond_price_lattice(11,1) = 1; 

  for col = [2:11]
    for row = [11:-1:1]
      if (col + row >= 12)
         
        if (row = 2)
          bond_price_u = q*bond_price_lattice(row, col-1)/(1+short_rate_lattice(row - 1, col - 1)/100);
        else
          bond_price_u = 0;
        endif  
        
        bond_price_lattice(row, col) = bond_price_u + bond_price_d;
      
      endif
    endfor 
  endfor

  for col = [2:11]
    bdt_model_zcb_prices(col - 1) = sum( bond_price_lattice(:, col) );
    bdt_model_spot_rates(col - 1) = 100 * (((1 / bdt_model_zcb_prices(col - 1)) ^ (1/(col -1)) ) - 1);
    squared_diffs(col - 1) = (bdt_model_spot_rates(col - 1) - spot_rate(col - 1)) ^ 2;
  endfor

  objective = sum(squared_diffs);
endfunction

I vet aquí la invocació a sqp

% main
clc

a = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3];

[x, obj, info, iter, nf, lambda] = sqp(a, @bdt_objective);

x'
obj
info
iter
nf
lambda

El resultat amb octave és, per cert, força més precís. I d’una eficiència similar a la de l’Excel.

albert@athena:~/MOOC/MOOC Financial Engineering and Risk Management - 1/Week5$ octave
GNU Octave, version 3.6.2
Copyright (C) 2012 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  For details, type `warranty'.

Octave was configured for "x86_64-pc-linux-gnu".

Additional information about Octave is available at http://www.octave.org.

Please contribute if you find this software useful.
For more information, visit http://www.octave.org/help-wanted.html

Read http://www.octave.org/bugs.html to learn how to submit bug reports.

For information about changes from previous versions, type `news'.

octave:1> ps5
ans =

   3.0000   3.1202   3.2326   3.3377   3.4357   3.5269   3.3095   3.3113   3.3111   3.3089

obj =  6.4126e-16
info =  101
iter =  45
nf =  46
lambda = [](0x1)
octave:2>