29 jul 2014

Análisis de la volatilidad en los mercados financieros usando Matlab


1.       Visión descriptiva de la volatilidad usando índices 

La volatilidad, junto al rendimiento y la liquidez forman variables básicas al momento los agentes tomar decisiones orientadas a optimizar sus decisiones. En los mercados financieros la volatilidad suele medirse a partir de las desviaciones o varianza del cambio porcentual del índice o precio (rendimientos), aunque esta medida asume varios retos que veremos brevemente adelante. En finanzas existen diversos tipos de riesgo: riesgo de mercado, de reinversión, operacionales, de crédito… por tanto la medición del riesgo requerirá definir con anterioridad que se desea medir y como. Además en el caso de los activos específicos cuentan con riesgos sistemáticos que se recogen dentro de un índice de mercado (no diversificable dentro del mismo mercado) y un riesgo no sistemático o propio de cada activo.

Ejemplo 1. Calculo de la volatilidad sobre el índice DowJone
En el siguiente ejemplo se utiliza la cotización de cierre diaria, del DowJone correspondientes a los meses abril-junio de 2014. Luego de cargar los datos se calcula la tasa de crecimiento o rendimiento del activo sobre el cual se calcula la varianza y desviación, que representan las volatilidad del indicice.[1]

indice = xlsread('DowJoneAbr14Jun14.xlsx')
plot(indice(:,1))
    title('Grafico 1. Evoulución del Dow Jone, Abril-Junio 2014')

format bank
  indice = xlsread('DowJoneAbr14Jun14.xlsx')

rendimiento = 100.0*(indice(2:end,1)./indice(1:end-1,1)-1.0)
  var(rendimiento)  %Varianza
  std(rendimiento)  %Desviación Estándar

Se pueden realizar algunas notas aclaratorias sobre el comando anterior, [xlsread] permite cargar la base de datos al sistema y [rendimiento] calcula el rendimiento del índice como el cociente entre el valor presente y el valor pasado, basado en dos vectores tomado sobre el índice el primer vector captura los valores actuales y va de la segunda observación a la última en la primera columna de la matriz índice [indice(2:end,1)], el segundo vector va desde la primera observación hasta la penúltima en la misma matriz [indice(1:end-1,1)] y representa los valores retardado. A este cociente le resta uno tal como se hace común mente.

El próximo comando (se anexa al final del documento), aunque un poco cutre, permite obtener una tabla con un resumen estadístico de los índices ('DJONE', 'SP500', 'IBMEX', 'NIKKE', 'IBEX35')

Tabla 1. Resumen estadístico valores de índice a nivel
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'Estadístico'           'DJONE'       'SP500'      'IBMEX'        'NIKKE'       'IBEX35'
---------------------------------------------------------------------------------------- 
    'Media'          [16923.31]    [1958.98]    [ 43069.15]    [15231.53]    [10916.05]
    'Varianza'       [14459.18]    [ 326.15]    [514025.52]    [20656.12]    [46200.55]
    'Dst.Est.'       [  120.25]    [  18.06]    [   716.96]    [  143.72]    [  214.94]
    'Coef.Var.'      [  140.74]    [ 108.47]    [    60.07]    [  105.98]    [   50.79]
    'Asimetria'      [   -0.07]    [  -0.59]    [    -0.14]    [   -0.33]    [   -0.56]
    'Curtosis'       [    2.03]    [   2.34]    [     2.47]    [    1.85]    [    2.09]
----------------------------------------------------------------------------------------

Es importante recordar que una mayor varianza no permite, por si sola, definir cuál de los índices ha sido más volátil. Sin Embargo el uso de desviaciones típicas como porcentaje de la media ['Coef.Var.'] si permite realizar este tipo de argumentos, por el hecho de expresar todos los valores en medidas comparables. Para obtener los valores en rendimientos, ahora no de un vector, sino de una matriz completa de datos, se modifica levemente el comando anterior y se puede obtener la misma tabla, pero sobre los rendimientos y anexando el test de normalidad de Jarque Bera ['Jarque-Bera'].[2]

rendimientos = 100.0*(indices(2:end,:)./indices(1:end-1,:)-1.0)

Observe que la única diferencia con el comando anterior de rendimientos es que se sustituye el 1 por los dos puntos (parte sombreada) esto le indica al programa que no solo realice el cálculo en la primera columna, sino que lo repita en todas. Otros lenguajes de programación como R, permiten hacer este tipo de operaciones. Una vez calculado el rendimiento de todos los índices se repite la tabla anterior ahora sobre el rendimiento.

Tabla 2. Resumen estadístico valores de índice en rendimiento
--------------------------------------------------------------------------
    'Estadístico'    'DJONE'    'SP500'    'IBMEX'    'NIKKE'    'IBEX35'
--------------------------------------------------------------------------
    'Media'          [ 0.07]    [ 0.09]    [ 0.20]    [ 0.05]    [ -0.08]
    'Varianza'       [ 0.19]    [ 0.17]    [ 0.39]    [ 0.42]    [  0.60]
    'Dst.Est.'       [ 0.44]    [ 0.41]    [ 0.63]    [ 0.65]    [  0.77]
    'Coef.Var.'      [ 0.15]    [ 0.23]    [ 0.32]    [ 0.07]    [ -0.10]
    'Asimetria'      [-0.73]    [-0.44]    [ 0.61]    [ 0.34]    [  0.04]
    'Curtosis'       [ 3.26]    [ 3.28]    [ 3.56]    [ 2.95]    [  3.10]
    'Jarque-Bera'    [ 0.00]    [ 0.00]    [ 0.00]    [ 0.00]    [  0.00]
--------------------------------------------------------------------------

Es importante resaltar que el cuadro anterior representa datos de un periodo de tiempo relativamente corto, por lo que la volatilidad puede estar afectada por una condición particular a ese periodo y no a condiciones estructurales de los mercados, por lo que para periodos suficientemente largos, las conclusiones pueden cambiar considerablemente. En ocasiones también es necesario obtener estimación de la volatilidad, por ventadas. Matlab permite obtener este usando un bucle que vaya rotando la ventana de acción y se ejecute por cada columna.

n = 15.0 %tamaño de la ventana;
for j=1:5
     for i = 1:((length(indices)-n))
       ventana_std(i,j) = std(indices(i:i+n,j));
      end
    end

ventana_std

El comando anterior[3] es un bucle anidado que corre el primer for y luego el segundo. En palabras el primer for [for j=1:5] le indica al programa que camine por columna y luego el segundo por observaciones [for i = 1:((length(indices)-n))] teniendo en cuenta la amplitud de la ventada. Es decir que en la primera secuencia de cálculo, j=1 indicando la columna e i estará caminando desde 1 hasta el máximo del vector menos la longitud de ventanas deseada (n=15), en este caso de 15 por lo que ira hasta end-15 observaciones. Ya dentro del bucle para la primera secuencia, la sentencia de cálculo de la volatilidad [ventana_std(i,j) = std(indices(i:i+n,j))] guarda dos elementos claves [ventana_std(i,j)] que permite que los datos se guarden en una matriz y no uno encima de otro y [indices(i:i+n,j)] que se va moviendo en el rango sobre el cual se calcula la volatilidad, en la primera observación será (1:1+15, 1) es decir en una ventana con las primeras 15 observaciones, luego (2:2+15, 1)… hasta que recorre todas las opciones de la primera columna, luego lo repetirá pero ira cambiando el 1 por 2,3 así hasta 5 indicando el cambio de columna.


En ocasiones es importante que estos datos sean los más comparables posibles, por ende es interesante comparar (en ventanas, como en el comando anterior) el coeficiente de variación, ya que presenta la desviación en proporción a la media, para obtenerlo solo habría que modificar la línea [4] sustituyéndola por la fórmula:

ventana_coefVar(i,j)= std(indices(i:i+(n-1),j))./mean(indices(i:i+(n-1),j));

Observe que la idea se repite, lo importante es especificar la forma en cómo se ira moviendo los rangos de datos. Siguiendo la idea de hacer las variables comparables es frecuente que estas se tipifiquen o se normalicen que no es más que restar la media a cada observación y dividir entre la desviación estándar. Este procedimiento se puede realizar en Matlab con relativa facilidad.

indices = xlsread('indicesJun14Jul14.xlsx',1)
format short
des.sta=std(indices)

for j=1:5
 for i = 1:length(indices)
        normalize(i,j) = indices(i,j) - mean(indices(:,j));
        normalize_indices(i,j) = normalize(i,j)/ des.sta(j);
     normalize_indices
  end
end

En caso que se desee llamar un vector de dentro de una matriz, no se necesita usar bucle y sustituyen la [i=:] y [j=nocol] es decir j será igual al número de columna que se desea extraer. En el caso de un vector la notación se simplifica considerablemente.

x=1:10;
xtipificada = x-mean(x)./std(x)


Este tipo de variables además de utilizarse para hacer las series comparables, es útil al momento de realizar análisis de correlaciones. Cuando se presentan los datos en función de la media se resuelve un problema importante en términos de comparabilidad aunque sigue presente el problema de la no estacionariedad de los índices. Por lo que es necesario trabajar con una transformación que si lo sea, por lo general la rentabilidad logarítmica. Para obtener la rentabilidad logarítmica de un vector o matriz se realizan los siguientes cálculos.

%En el caso de un vector
DJONE = indices(:,1);
  rlogDJONE = log(DJONE(2:end)./DJONE(1:end-1))*100.0

%En el caso de una matriz de datos
 Rlog_indices = log(indices(2:end,:)./indices(1:end-1,:))*100.0

La volatilidad se comúnmente utilizada para para determinar estrategias de cobertura, por tanto definirla mal daría lugar a estrategias inadecuadas, además del uso extendido de la matriz de varianza covarianzas en la gestión de riesgo. Por lo que es importante tener buena medida de la volatilidad y considerar las limitaciones de la estimación de la volatilidad muestra. Entre esta se encuentra: el sesgo al estimar las desviaciones, la posible presencia de tendencias deterministas o estocásticas y cualquier cambio de nivel en las variables hace perder representatividad a las medidas de tendencia centrales.
Por ultimo dos comandos que suelen ser útiles en una gran cantidad de aplicaciones financieras son la matriz de correlaciones y la de varianza covarianzas. Ambas son cuadradas y simétricas. En el caso de la primera se obtiene con el comando [cov] y arroja una matriz donde la diagonal principal muestra las varianzas de cada serie y los elementos fuera de la diagonal principal la covarianza entre las series. En tanto que la matriz de correlaciones se llama usando [corr] esta muestra la correlaciones entre las series, siendo la diagonal principal igual a 1 porque representa la relación lineal de la serie con ella misma.

cov(indices)

    0.1446    0.0186    0.7690    0.1050   -0.1142
    0.0186    0.0033    0.1108    0.0205   -0.0102
    0.7690    0.1108    5.1403    0.5547   -0.7333
    0.1050    0.0205    0.5547    0.2066   -0.0347
   -0.1142   -0.0102   -0.7333   -0.0347    0.4620

corr(indices)

    1.0000    0.8568    0.8920    0.6077   -0.4418
    0.8568    1.0000    0.8555    0.7897   -0.2639
    0.8920    0.8555    1.0000    0.5383   -0.4759
    0.6077    0.7897    0.5383    1.0000   -0.1125
   -0.4418   -0.2639   -0.4759   -0.1125    1.0000

2.      Medición de la volatilidad de un activo financiero

Se puede considerar el cuadrado de la rentabilidad diaria para utilizarlo como predicción de la volatilidad diaria. En Matlab este se haría de la forma siguiente, primero se estima el rendimiento como la diferencia logarítmica y luego cada uno de estos valores [.^] se eleva al cuadrado.

DJONE = indices(:,1);
  R_Djone = log(DJONE(2:end))-log(DJONE(1:end-1));  
  Vol_Djone = R_Djone.^2 ;

plot(t, Vol_Djone);
  
1.1. Medidas de Parkinson y Garman-Klass

Frente a las deficiencias encontradas en la volatilidad como medida de riesgos se han incorporado otras medidas que pueden utilizarse dependiendo del tipo de información que se maneje. Cuando se dispone del precio de todas las transacciones diarias de un activo, la formula difiere a cuando no se tratan de transacciones diaria o cuando se dispone de pocas observaciones. En nuestro caso disponemos de transacciones diarias (pero la usaremos a modo de ejemplo como si se tratasen de transacciones hechas en un mismo días, por tanto la volatilidad anual se obtiene de la siguiente forma:

En la formula anterior la primera raíz sirve para anualizar la volatilidad y su valor puede variar según la transformación en el tiempo que se necesite (semanal=52), la segunda raíz permite calcular la volatilidad diaria usando las variaciones del rendimiento de las operaciones registradas durante ese día. Además existen otras propuestas para cuando se disponen de otro tipo de información. Siendo r el rendimiento. En el caso de la primera raíz, cuando se trate de una volatilidad diferente a diaria, se puede modificar para obtener la volatilidad anual. Ejemplo 12 en caso que sea mensual, r es igual a la rentabilidad diaria y T las observaciones.

Ejemplo.
Disponiendo de los datos mensuales de cotización diaria (aquí se entienden como si fueran todas transacciones hecha en un mismo día), en fecha de 07/05/2007 al 28/02/2014[4]. Una vez cargada la base solo habría que estimar el rendimiento y luego colocar la formula mostrada anteriormente. (Esta fórmula se puede ir modificando según las distintas propuestas)

%Cálculo de la rentabilidad
   r = log(azucar(2:end)./azucar(1:end-1))*100.0

%Cálculo de la volatilidad
   Vol_azucar = sqrt(12)*sqrt((1./length(r)-1).*sum((r-mean(r)).^2))

%siendo T = length(r)

1.1. Volatilidad implícita Vs. Volatilidad histórica

La volatilidad implícita se calcula a partir del precio de los activos, utilizando la ecuación de Black-Scholes que se presenta a continuación. Siendo c el precio de mercado de una opción de compra  europea y estando todos los parámetros conocidos, salvo sigma, lo que se intenta es buscar el valor de sigma que iguala el precio de mercado con el precio teórico. Este valor de sigma es lo que se conoce como volatilidad implícita. Es decir, se puede determinar la volatilidad implícita a partir del precio de la opción en ausencia de arbitraje, por lo que depende del deseo de compra y venta.

De forma llana se tiene una igualdad entre la formula anterior y el precio de mercado, solo falta el valor de sigma que hace que ambos precios se igualen.

Ejemplo. Calculo de la volatilidad implícita de una Call Digital
En el siguiente ejemplo se muestra como calcular la volatilidad implícita usando la fórmula de BS de una opción europea binaria. (Este mismo procedimiento se puede aplicar para casi cualquier opción, modificando los parámetros de la fórmula de cálculo, el payoff de la opción). En este caso se crean tres funciones que se expone a continuación.

Paso 1. La idea de las dos primeras funciones es calcular el precio de mercado de una opción que pague a en caso que el precio del subyacente sea mayor al precio de ejercicio (ST>K) y nada en caso contrario (payoff). (Este payoff se podría medicar según el tipo de opción sobre la cual se esté trabajando)

function f = integrandoCallDig(x,S0,K,r,T,sigma,a)
 %Calcula el precio del subyacente en el momento T
     ST = S0*exp((r-0.5*sigma^2)*T + sigma*sqrt(T)*x);
    
 %Calcula el Payoff (funcion de pago) de una Call Digital
        for i=1:length(ST)
           if ST(i)>=K 
              payoff(i) = a;
           else
              payoff(i) = 0;
          end
         end
   
   f = normpdf(x).*payoff;
end

function precio = precioCallDig(S0,K,r,T,sigma,a)
%precioCallDig: Precio de una call europea Digital
  precio = exp(-r*T)*quadl(@(x)integrandoCallDig(x,S0,K,r,T,sigma,a),-10,10);
end

En el caso que se estén usando simulaciones por Montecarlo se utiliza la cuadratura por Montecarlo para sustituir la integral de esta última formula. (En otro documento se explicara más adelante como se hace).

Paso 2. En el segundo paso se obtiene la estimación de la volatilidad implícita mediante un bucle while, la idea es ir modificando el valor de la volatilidad, en la función que calcula el precio teórico (segunda función), hasta que la diferencia entre el precio teórico y el precio del mercado sea igual a cero.

function f = Vol_Imp_CallDig(market,S0,K,r,T,sigma,a)
 %Vol_Imp_CallDig: Volatilidad implicita de una call europea Digital

 %Se especifica el formato de salida y se especifican los valores de
 %partidas de la diferencia entre el precio teórico y el precio de mercado.

  format long
    pBS = precioCallDig(S0,K,r,T,sigma,a);
    V   = pBS-market;

%Se crea un bucle (while)que modifica la volatilidad hasta que el precio de
%la Call obtenido por Black-Scholme, sea igual al precio de mercado

while 1
    if V<=0
        disp('A esta volatilidad los precios se igualan')
        break
    end
      sigma=sigma+0.0001;
        V = precioCallDig(S0,K,r,T,sigma,a)-market;
      f=sigma;
end

Una vez creadas y guardadas las funciones en el path o directorio de Matlab en el que se esté trabajando, se utiliza esta última función para obtener la volatilidad implícita:

Vol_Imp_CallDig(7.81,100,90,0.1,2,0.001,10)

ans =

    0.122000000000002

Ahora bien para graficar la volatilidad, el método recursivo (While) usado anteriormente, solamente funciona en caso que el vencimiento o los parámetros ofrezcan un precio teórico superior al precio de mercado. Por tanto con el ejemplo anterior no se puede utilizar vencimientos superiores a 3. Por tal motivo se modifica (de forma un tanto cutre) el comando, incluyendo un if anidando (comandos condicionales) que cambie la forma de igualar ambas expresiones dependiendo de si el valor teórico está por encima o por debajo del de mercado.

Una idea para graficar la volatilidad implícita a distintos vencimientos es crear un bucle que vaya cambiando los valores de T en la fórmula de vencimiento, así se puede obtener la volatilidad implícita vigente en cada vencimiento.

S0 = 100; K = 90; r = 0.1; T = 2; sigma = 0.001; a=10; market=7.81;

   for T = 1:10 %1:2
     Vol_implicita(T) = Vol_Imp_CallDig(market,S0,K,r,T,sigma,a);
  end

plot(T, Vol_implicita)



[1] Los datos usados le corresponden a las 35 observaciones siguientes al 2 de junio de 2014 de cada datos, cotizaciones de cierre.

[2] Se hace que la matriz indice se llame rendimiento para poder utilizar el mismo comando anterior [indices=rendimientos]

[3] Los bucle anidados para completar matrices se explica en un documento anterior, titulado matrices y tablas en Matlab.

22 jul 2014

El Proceso Browniano con Matlab

1. Proceso de Wiener o movimiento Browniano en Matlab

Una definición teórica del proceso de puede encontrar en: http://pendientedemigracion.ucm.es/info/jmas/mon/28.pdf

El movimiento representa por:

Que se genera y gráfica por el código:
-----------------------------------------------------------------------------------------------------
%Generando procesos de Wiener
format short
   B0=0; N=200; M=100; mu=0.1; sigma=0.4; T=2; DeltaT=T/N;

X = randn(M,N);    %X~N(0,1)
Bt = cumsum([B0*ones(M,1) (mu*DeltaT*sigma*sqrt(DeltaT)*X)],2);
 
  figure(1); plot(0:DeltaT:T, Bt', '-' ),
      title ('Representación de 200 trayectorias de Wiener');           
      xlabel('Tiempo');


2.   Modelo de Black-Scholes, movimiento Browniano geométrico en Matlab

Una definición teórica del proceso de puede encontrar en: 
http://es.wikipedia.org/wiki/Movimiento_browniano

El movimiento se representa por:


Que se genera y gráfica por el código:
------------------------------------------------------------------------------------------------------
%Generando procesos Browniano Geométrico
format short
S0 = 100.0; N = 1e2; M = 1e3; mu = 0.5; sigma = 0.4; T = 2; dT = T/N;

   X = randn(M,N);
  St = cumprod([S0*ones(M,1) exp((mu-0.5*sigma^2)*dT+sigma*sqrt(dT)*X)],2);

   figure(2); plot(0:dT:T, St')
hold on; plot(0:dT:T, mean(St),'k','linewidth',4); hold off
    title ('Representación de 200 trayectorias Browniano Geométrico');           
    xlabel('Tiempo');








19 jul 2014

Matrices y tablas en Matlab


3.1. Notas introductorias 

En esta parte se tratan matrices, algunas tablas y se ilustra cómo usar bucles para completar matrices. Aunque existe abundante literatura en la red. Las formas básicas de crear columna y matrices en Matlab es:

%Vector fila
>> A=[1 2 3]

A =

     1     2     3

%Vector columna
>> B=[1; 2; 3]

B =

     1
     2
     3

%Es equivalente
>> B = [1 2 3]'
>> B = [1:3]'

Estos dos últimos comando son equivalentes, en el primero transpone un vector fila en columna ( ' ), en el segundo se ilustra como usar operaciones vectorizadas en Matlab, lo que evita el uso de bucle cuyo uso para construir matrices se presenta en el apartado 3.6.

%Matriz bidimensional
C=[1 2; 2 4; 3 6]

C =

     1     2
     2     4
     3     6

Esta matriz bidimensional, seria un problema si se tratara de una matriz de mayor dimensión, sin embargo Matlab permite obtenerla de forma vectoriada. la matriz C, se puede obtener usando:

>> C = [1:3; 2:2:6]'

C =

     1     2
     2     4
     3     6

La primera parte del comando anterior ordena una secuencia de números que vaya desde el 1 hasta el 3, y una segunda columna que llama, también de forma vectorizada, a crear una columna, que vaya del 2 al 6, en pasos de 2 en 2. Un comando que realiza una operación similar a esta última citada es linspace, este permite dividir intervalos según fuese requerido. Por ejemplo, dividir en 5 un intervalo que va desde 0 hasta 10.  Es decir obtener un vector con x valores igualmente espaciados.

>> D = linspace(0,10,5)

D =

         0    2.5000    5.0000    7.5000   10.0000

Nótese que el uso de los dos puntos (:) tiene un uso extendido en Matlab. Se puede acceder a la primera columna de una matriz especificando (:,1), accedería a la columna 1.

>> C(:,1)

ans =

     1
     2
     3

Si se usa C(:,2:end) recorrerá todas las filas, desde la segunda hasta la última columna. El siguiente comando permite representar conjuntamente, datos de números aleatorios que se generaron de forma independientes. Matlab también permite representar conjuntamente datos que han sido generado de forma independiente. El siguiente comando permite representar conjuntamente, datos de números aleatorios que se generaron de forma independientes. (el ; se coloca para poder colocar varias instrucciones en una misma linea) 

>> y = [rand(1,10)]';, x=[rand(1,10)]';, z=[y,x]

z =

    0.7513    0.8407
    0.2551    0.2543
    0.5060    0.8143
    0.6991    0.2435
    0.8909    0.9293
    0.9593    0.3500
    0.5472    0.1966
    0.1386    0.2511
    0.1493    0.6160
    0.2575    0.4733 

En resumen y dado a la abundante literatura respecto al tema que existe en la red, en la próxima tabla se resume un conjunto de comandos básicos de operaciones con matrices:

Construir matrices

Forma de una matriz
matriz(fila, columna)
dimensiones de la matriz
[nfilasM, ncolsM]=size(A)
    Crear una matriz
A = [16 3 2; 1 11 8; 6 7 1; 4 1 1]
    matriz a ceros
A = zeros (3, 5)
    matriz a unos
A = ones (3, 5)
    matriz identidad
eye(3)  
    Matrices x distantes
linspace(xmin, xmax, n)
x=xmin:stepsize:xmax


Llamar o editar un elemento
A(1,4)
  Última fila
sum(A(:,end))
  Borrar la segunda columna
X(:,2) = []


builds a table
x = (1:0.1:2)’;
logs = [x log10(x)]


Matrices esperciales

   Transpuesta
A’
   matriz a ceros
A = zeros (3, 5)
   matriz a unos
A = ones (3, 5)
   matriz identidad
eye(3)  
   dimensiones de la matriz
[nfilasM, ncolsM]=size(A)
   Tomar la diagonal
diag(A)
   The determinant
d = det(A)
   Inversa
X = inv(A)
   Triangulares
tril(A) triu(A)
   Autovectores y Autovalores
[V, D] = eig(A)
   Cholesky
chol(A)
   Valores singulares
svd(A)


Operaciones con matices

Sumar

   Sumar filas
sum(A’)’
   Sumar la diagonal
sum(diag(A))
   Crear matriz simétrica
A + A’
multiplicación

   Multiplicar
A’*A
Exponencial matricial
expm(A)


Estadísticos

Estadísticos por columna
mu = mean(D);
sigma = std(D)
Nota. En cuanto a los estadísticos de matices, en la parte de Estadística descriptiva se explicó su uso para trabajar por fila o columna.

Trabajar con matrices abre muchas posibilidades, por ejemplo podemos generar una matriz con los números pares entre 1 y 10 y sus respectivos logaritmos. Este ultimo comando es especialmente útil para cuando se desea calcular factoriales dobles de forma vectorizada con Matlab.

>> num_pares = (2:2:10)'
>> Log=[num_pares log(num_pares)]

Log =

    2.0000    0.6931
    4.0000    1.3863
    6.0000    1.7918
    8.0000    2.0794
   10.0000    2.3026

3.2. Poner texto y números en un cuadro

Esta tarea no es trivial en Matlab. No es posible insertar textos en las matrices de Matlab, si se llaman directamente. Se deben usar comillas simples y llamar mediante llaves (no corchetes).

>> nom={'nerys', 'salami', 'platano'}
>> class(nom)

Este último comando permite utilizar   “cell”, comando que permite guardar diversos tipos de datos. Cuyos objetos y dimensiones se acceden mediante {}. En el siguiente ejemplo se muestra como especificar una tabla con los nombres y la media de un conjunto de datos.

% Ejemplo: Clasificacioens promedio por notas obtenidas
%Especificas clases y datos
nota  = [58 67 80 96]
notas = {'Rep','Extra','Apro','Sobre'}

%Completando tablas
  cellnotas=cell(4,2);
     cellnotas (:,1)=notas; % los : indica es una operación vectorizada

       %se pueden además completar celdas directamente
       Rep   = nota<70; Extra = nota<70 & nota>60;  
       Apro  = nota>70; Sobre = nota>85;

    cellnotas(1,2)= {mean(Rep)};,   cellnotas(2,2)= {mean(Extra)};
    cellnotas(3,2)= {mean(Apro)};,  cellnotas(4,2)= {mean(Sobre)};
cellnotas

cellnotas =

    'Rep'      [0.5000]
    'Extra'    [0.2500]
    'Apro'     [0.5000]
    'Sobre'    [0.2500]

3.3. Clasificaciones

En ocasiones es necesario crear tablas donde se incluyan además el nombre de la clase con que se esté trabajando. Por ejemplo se pueden hacer tablas de frecuencia, incluyendo clasificaciones específicas.

%Tabla de frecuencia con etiquetas de clases
>> nota  = [75 85 92 15 65 82 72 45 94 85 86 70 68 56 57 63 72 14]
>> notao = ordinal(nota,{'Rep','Extra','Apro','Sobre'},{},[0,59,69.9,85,100])

>> tabulate(notao)

  Value    Count   Percent
    Rep        5     27.78%
  Extra        3     16.67%
   Apro        5     27.78%
  Sobre        5     27.78%

3.4. Tablas y matrices con operadores lógicos

El uso de operadores lógicos se supone conocido. Utilizar operadores lógicos para obtener información de las tablas y vectores calculados en Matlab, se presenta a continuación. 

%Calcular la media de estudiantes susmendidos
nota  = [75 85 92 15 65 82 72 45 94 85 86 70 68 56 57 63 72 14]
      suspenso = nota<70         %indica 1/0. 1 si suspendieron               
     nota_susp =nota(suspenso)  %muestra la nota de los que suspendi…
      mean(nota_susp)            %media de notas de los que suspedieron

El promedio de estudiante suspendidos, se puede obtener, aprovechando que la matriz [suspenso] es dicotómica o binaria, llamando el promedio de dicha matriz o calculando directamente la media, el próximo cuadro presenta respectivamente, dichos comandos.

>> mean(suspenso)
>> sum(suspenso)/length(nota)

Existen en Matlab un conjunto de operadores lógicos que permiten aplicar pruebas lógicas a matrices de datos. 

Comandos
Descripción
any(nota>85)
%arroja un 1 si en la matriz nota hay por lo menos un número mayor a 85 (cumpla alguna condición)
any(C(:,1)>=3)
%Aplica la misma condición sobre una fila
all(nota>85)
%revisa si todos los elementos cumplen una condición
find(nota>85)
Busca valores y posiciones que cumplan determinadas condiciones

Ejemplo: usar la función find para sustituir los elementos de una matriz que cumplen determinadas condiciones.

N=magic(4)

N =

    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1

%buscar valores mayores a 10
>>    J = find(N>=10)
>> N(J) = 10*ones(size(J))

N =

    10     2     3    10
     5    10    10     8
     9     7     6    10
     4    10    10     1

         3.5. Tablas de contingencia y etiquetas filas y columnas

Teniendo dos vectores de datos, X referido al sexo 1==hombre 0==mujer y una segunda variable fuma que muestra un 1 en el caso que la persona fume. El comando usado agrega el chi2 y el pvalor, usando la hipótesis nula de independencia entre variables.

>> sexo = round(rand(1,10))'
>> fuma = round(rand(1,10))'  
>> [table,chi2,p] = crosstab(sexo, fuma)

table =

     2     4
     2     2


chi2 =

    0.2778


p =

    0.5982

Hasta ahora hemos obtenido todas las tablas y matrices sin las debidas etiquetas, que permitan identificar su contenido. La tabla de contingencia anterior se podría modificar de forma que aparezcan las debidas etiquetas.

%ahora modificamos el comando para obtener las debidas etiquetas
table = array2table(table,              ...
      'VariableNames',{'Hombre','Mujer'},...
      'RowNames',{'Fumador','No_fuma'});

disp('Tabla cruzada')
 disp(table)

Tabla cruzada
                 Hombre       Mujer         
                --------    ---------   
    Fumador            2            4    
    No_fuma            2            2   

3.6. Usar bucles para completar matrices de datos

Una forma adicional de construir matrices se encuentra en los bucles. Por ejemplo, la Matriz B, obtenida en el primer apartado, se podría obtener alternativamente usando:

>>  for I = 1:3,
      n(i) = i;
     end

>>  n =

     1     2     3

Esta forma de construer matrices se vuelve especialmente útil cuando necesitamos completar matrices bidimensionales, siguiendo el criterio de completar ya sea por fila o por columna. 

>> %Recorrer una matriz por fila
for i=1: nfilasA                 % para cada fila
     for j=1: ncolsA
              sentencias
       end
 end

% recorrer la matriz por columnas
for j=1: ncolsA                  % para cada columna
      for i=1: nfilasA
              sentencias
        end
end

en el primer caso se tiene:

 Si aun no se ve se puede aclarar viéndolo arbitraria donde cada numero muestra en que paso el programa trabajo en esa celda. en los próximo se muestra como el programa va trabajando según el criterio que se va indicando. En el primero ejecuta la sentencia en la primera fila y trabaja todos sus elementos columna por columna y en la segunda trabaja todos los elementos de la columna fila por fila.

>> %Recorrer una matriz por fila
n =

     1     2
     3     4
     5     6

recorrer la matriz por columnas
n =

     1     4
     2     5
     3     6

Para ver una aplicación un poco mas avanzada de estos comando, en el siguiente comando se utiliza el bucle anidado for para simular un proceso autoregresivo (AR(1)) de primer orden, rellenando una matriz por fila. La clave esta en el orden en que se colocan los for, pues estos indican la forma en como deberá ir trabajando el programa.

M = 2; T = 12; X0 = 10; phi0 = 0.0; phi1 = 0.95; sigma = 0.25;
X= zeros(M,T);

%Se completa la matriz anterior con un bucle que completa, por fila, el
%numero de trayectorias especificadas en la función.
   for i=1:M                             
       for j=2:T                             
           u(i,j)=sigma*rand();
          X(i,1)=phi0+phi1*X0+rand()*sigma;
          X(i,j)= phi0+phi1*X(i,j-1)+u(i,j);
       end
   end
  

REFERENCIAS BIBLIOGRÁFICAS

Lopez, Beatriz (2009). Estadística. Disponible en: http://eio.usc.es/pub/pateiro/files/iq0809pateiro.pdf

Creando variables por grupos en dplyr (group_by + mutate)

  Simulemos una base de hogares, donde se identifica el hogar, el sexo (1 mujer) y provincia y edad para cada miembro.   # Definir la lista ...