Introducción aplicada del uso de Matrices en R: 9 aplicaciones en economía

Borrador...

De manera puntual, la siguiente entrada toma elementos relacionados con las matrices, de los libros: Análisis econométrico de William Greene; Métodos Fundamentales de Economía Matemática de Chiang y Wainwright, 2006; e Introducción a la econometría de Wooldridge (4a. edición). Lo anterior, para intentar explicar algunas definiciones elementales de matrices, acompañada de alguno de sus usos más comunes en economía, como es el caso de resolución de sistemas de ecuaciones a partir de una matriz inversa o la determinación de componentes principales a partir de la matriz de varianzas-covarianzas de la serie de datos normalizada.

Estos ejemplos se realizan en el programa R, intentando motivar la importancia del uso de matrices en economía en aspecto como los modelos de transición de estado, optimización de funciones o resolución de problemas de cartera en el campo de la economía financiera, donde (en el caso de las finanzas), además de obtener los ponderadores de la cartera de mínima varianza sin necesidad de resolver el problema de optimización, puede obtenerse la varianza de una cartera determinada, solo teniendo los ponderadores de los activos que confirman la cartera y su matriz de varianza covarianza y utilizando el producto de matrices, lo que ahorra tener que utilizar la expresión algebraica del cálculo de la varianza de activos (que puede ser bastante larga en cartera de muchos activos).

Vectores y matrices

Una matriz puede verse como un conjunto de vectores columnas (Greene, p.18). Los números dentro del vector son llamados componentes o coordenadas del vector. El componente vi es llamado el i-ésimo componente del vector. Es decir, qeu un vector con n componentes, es llamado n-vector. Las matrices se ordenan en filas y columnas, cuyas longitudes, definen las dimensiones.

En R, una alternativa para crear matrices es utilizar el comando matrix, cuyos argumentos fundamentales son:

matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)

# vector renglón (fila)
vector1<-c(1,8,4,-1)
vector1

## [1]  1  8  4 -1

matriz1 <- matrix(vector1, 2, 2)
matriz1

##      [,1] [,2]
## [1,]    1    4
## [2,]    8   -1

La indexación de matrices en R, se especifica indicando las finas y columnas a las cuales se desea acceder, y aplicando los criterios de indexación tradicionales:

#A11
matriz1[1,1]

## [1] 1

#Toda una columna
matriz1[ ,1]

## [1] 1 8

# Toda una fila
matriz1[2, ]

## [1]  8 -1

Matriz simétrica

Una Matriz Simétrica es aquella matriz cuadrada que es igual a su traspuesta (Chiang y Wainwright, 2006, p.72). Por lo que, el elemento A[i,j]=A[j,i], siendo lo anterior valido para todos los argumentos de la matriz, por lo que, la simetría se produce alrededor de los valores de la diagonal principal de la matriz.

matriz2<- matrix(c(1,4,3,
                   4,2,5,
                   3,5,6), 3, byrow = F)

matriz2

##      [,1] [,2] [,3]
## [1,]    1    4    3
## [2,]    4    2    5
## [3,]    3    5    6

isSymmetric(matriz2)

## [1] TRUE

Matriz diagonal

Una matriz cuadrada es una matriz diagonal si todos sus elementos fuera de la diagonal son cero (Wooldrige, 2009, p.789). Observe que la función diag permite obtener una matriz diagonal de un determinado vector que se le introduce como argumento.

matriz3 <- diag(c(vector1,4))
matriz3

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    8    0    0    0
## [3,]    0    0    4    0    0
## [4,]    0    0    0   -1    0
## [5,]    0    0    0    0    4

Tenga pendiente que al utilizar la función diag dándole una matriz cuadrada como argumento, esta regresa un vector con los elementos de la diagonal principal de la matriz indicada.

diag(matriz3)

## [1]  1  8  4 -1  4

Matriz escalar

Una matriz escalar es una matriz diagonal en la que los elementos de la diagonal principal son iguales. 

matriz4 <- diag(rep(4,5))
matriz4

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    4    0    0    0    0
## [2,]    0    4    0    0    0
## [3,]    0    0    4    0    0
## [4,]    0    0    0    4    0
## [5,]    0    0    0    0    4

Matriz identidad

La matriz identidad se define como una matriz cuadrada (tiene el mismo número de filas que de columnas) con unos en su diagonal principal y ceros en cualquier otra parte (Chiang y Wainwright, 2006, p.70). Es decir, es una matriz escalar con 1 en su diagonal principal. Esta desempeña la función del número 1, de forma tal que: AI=IA=A.

matriz5 <- diag(rep(1,5))
matriz5

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    1    0    0    0
## [3,]    0    0    1    0    0
## [4,]    0    0    0    1    0
## [5,]    0    0    0    0    1

La traza de la matriz

Sea una matriz cuadrada M de orden n, la traza (tr(M)) es la suma de los elementos de la matriz (definida sólo para matrices cuadradas).

sum(diag(matriz5))

## [1] 5

Matriz triangular

Una matriz triangular es una matriz cuadrada cuyos elementos arriba o debajo de su diagonal principal son todos cero.

matriz6 <- diag(c(1,5,4,3,5))
matriz6[lower.tri(matriz6)] <- (1:10)
matriz6

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    1    5    0    0    0
## [3,]    2    5    4    0    0
## [4,]    3    6    8    3    0
## [5,]    4    7    9   10    5

lower.tri(matriz6)

##       [,1]  [,2]  [,3]  [,4]  [,5]
## [1,] FALSE FALSE FALSE FALSE FALSE
## [2,]  TRUE FALSE FALSE FALSE FALSE
## [3,]  TRUE  TRUE FALSE FALSE FALSE
## [4,]  TRUE  TRUE  TRUE FALSE FALSE
## [5,]  TRUE  TRUE  TRUE  TRUE FALSE

En el caso del uso de modelos de vectores autoregresivos, en la identificación recursiva de los coeficientes, el uso de este tipo de matrices adquiere un protagonismo especial para la reconocida descomposición de Cholesky.

Transpuesta

Es una nueva matriz que se obtiene, donde todas las columnas representan las filas de la matriz original, o viceversa.

t(matriz6)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    1    2    3    4
## [2,]    0    5    5    6    7
## [3,]    0    0    4    8    9
## [4,]    0    0    0    3   10
## [5,]    0    0    0    0    5

Note que ahora tenemos las herramientas para crear una matriz simétrica en R, utilizando la transpuesta y las matrices triangulares.

matriz6<-t(matriz6)
matriz6[lower.tri(matriz6)] <- (1:10)
matriz6

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    1    2    3    4
## [2,]    1    5    5    6    7
## [3,]    2    5    4    8    9
## [4,]    3    6    8    3   10
## [5,]    4    7    9   10    5

isSymmetric(matriz6)

## [1] TRUE

Propiedades:

(A’)’=A

(A+B)‘=A’+B’

(A*B)‘=B’A’

Matriz nula

Una matriz nula, o matriz cero, denotada por 0, juega el papel del número cero.

matriz5 <- diag(rep(0,5))
matriz5

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    0
## [2,]    0    0    0    0    0
## [3,]    0    0    0    0    0
## [4,]    0    0    0    0    0
## [5,]    0    0    0    0    0

Operaciones con matrices

Suma y resta

Note que suma los elementos homólogos colocados en cada matriz, por lo que, se requiere que las matrices tengan las mismas dimensiones:

A <- matrix(1:4, 2); B <- matrix(4:1, 2)
A; B

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

##      [,1] [,2]
## [1,]    4    2
## [2,]    3    1

A+B

##      [,1] [,2]
## [1,]    5    5
## [2,]    5    5

A-B

##      [,1] [,2]
## [1,]   -3    1
## [2,]   -1    3

Propiedades:

#Comutativa
all(A+B == B+A)

## [1] TRUE

all(t(A+B)==t(A)+t(B))

## [1] TRUE

Producto de matrices

Producto interno: dos matrices

a<-matrix(c(1,3,4), 3, 1); b<-matrix(c(3,8,2), 3, 1)
a; b

##      [,1]
## [1,]    1
## [2,]    3
## [3,]    4

##      [,1]
## [1,]    3
## [2,]    8
## [3,]    2

# Multiplicación elemento a elemento
sum(a*b)

## [1] 35

t(a)%*%b

##      [,1]

## [1,]   35

a’b = b’a

all(t(a)%*%b == t(b)%*%a)

## [1] TRUE

La multiplicación elemento a elemento (representada en R por el símbolo aritmético de multiplicación *), solo coincide con el producto interno de matrices en el caso de vectores. El operador para multiplicar matrices en R es %*%. En el caso de matrices, el número de columna de la primera, debe ser igual al número de fila de la segunda (conformabilidad para la multiplicación). Note que esta multiplicación suma el producto de los elementos correspondientes a la primera fila de la matriz A, vs. los argumentos de la primera columna de la matriz B.

# Greene, p.8
A <- matrix(c(1,4,3,5,2,1), 3, byrow = T)
B <- matrix(c(2,1,0,4,6,5), 3, byrow = F)
A; B

##      [,1] [,2]
## [1,]    1    4
## [2,]    3    5
## [3,]    2    1

##      [,1] [,2]
## [1,]    2    4
## [2,]    1    6
## [3,]    0    5

C<-t(A)%*%B
C

##      [,1] [,2]
## [1,]    5   32
## [2,]   13   51

Verifique que las siguientes funciones no corresponden a multiplicación interna. La segunda multiplica cada elemento de una matriz, con el matriz correspondiente de la segunda.

sum(A*B)

## [1] 56

A*B

##      [,1] [,2]
## [1,]    2   16
## [2,]    3   30
## [3,]    0    5

Transpuesta del producto

t(A*B)

##      [,1] [,2] [,3]
## [1,]    2    3    0
## [2,]   16   30    5

t(B)*t(A)

##      [,1] [,2] [,3]
## [1,]    2    3    0
## [2,]   16   30    5

all(t(A*B)==t(B)*t(A))

## [1] TRUE

Producto interno: matrices por escalar

Multiplica cada elemento por el escalar.

B*2

##      [,1] [,2]
## [1,]    4    8
## [2,]    2   12
## [3,]    0   10

Producto interno: matrices por vectores

B%*%c(1,2)

##      [,1]
## [1,]   10
## [2,]   13
## [3,]   10

Se lo anterior se deriva:

D*I = D

D <- matrix(c(1,2,3,4), 2, byrow = T)
D

##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4

I <- diag(rep(1,ncol(D)))

D%*%I

##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4

all(D%*%I == D)

## [1] TRUE

D*0 = 0

D%*%c(0,0)

##      [,1]
## [1,]    0
## [2,]    0

Aplicación 1. Suma de elementos de un vector

Verificar que una aplicación común del uso del producto de matrices y la matriz identidad, es obtener resumir la suma de un vector. Por ejemplo: t(I)%*%x, es equivalente a obtener la suma de los argumentos del vector x.

x<-c(1,6,4,3,1)

sum(x)

## [1] 15

I<-rep(1,length(x))
I

## [1] 1 1 1 1 1

t(x)%*%I

##      [,1]
## [1,]   15

t(I)%*%x

##      [,1]
## [1,]   15

En el caso de una constante a, su suma n veces es igual a n*a, se obtiene matricialmente mediante I’(aI). Por ejemplo, si queremos sumar a=4, 5 veces, utilizamos un vector de unos de longitud 5.

a<-4
t(I)%*%(a*I)

##      [,1]
## [1,]   20

sum(rep(a,5))

## [1] 20

Ahora, dado una constante y un vector x, la suma del producto de los valores del vector, multiplicado por la constante, a*sum(x), se obtiene mediante:

a*sum(x)

## [1] 60

a*t(I)%*%x

##      [,1]
## [1,]   60

Por tanto, si a=1/n, se obtiene la media del vector:

a<-1/length(I)

a*t(I)%*%x

##      [,1]
## [1,]    3

mean(x)

## [1] 3

En el caso de la suma de cuadrados de un vector, es necesario realizar la multiplicación interna de x’x.

sum(x^2)

## [1] 63

t(x)%*%x

##      [,1]
## [1,]   63

Producto externo

Finalmente, observe que x’x es diferente a xx’.

x%*%t(x)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    6    4    3    1
## [2,]    6   36   24   18    6
## [3,]    4   24   16   12    4
## [4,]    3   18   12    9    3
## [5,]    1    6    4    3    1

A*B’

Esta última operación se conoce como producto externo, ahora, cada fila de la matriz resultante corresponde al producto de los elementos de cada fila en la primera matriz, multiplicada por cada elemento de la columna correspondiente en la segunda matriz:

matriz8<-c(1,2,3)
matriz8%*%t(matriz8)

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    2    4    6
## [3,]    3    6    9

Aplicaciones 2. Gestión de carteras

En el siguiente ejemplo, tomado Tomado de Marín y Rubio (2010, p.203), se necesita estimar la varianza de una cartera de n activos (en el ejemplo siguiente n=3), donde se invierte w1=40%; w2=35%; w3=25%. Esta verianza se puede obtener realizando la multiplicación interna de w’(VarCov)w.

at<- c(12.75, 14.37, -15,

       14.37, 16.69, -16.5,

       -15, -16.5, 18)

MVar_Cov <-matrix(dat, 3,3)

MVar_Cov

##        [,1]   [,2]  [,3]

## [1,]  12.75  14.37 -15.0

## [2,]  14.37  16.69 -16.5

## [3,] -15.00 -16.50  18.0

w<-c(0.4,0.35,0.25)

Var_cartera<-t(w)%*%MVar_Cov%*%w

Var_cartera

##          [,1]

## [1,] 3.345625

Descomposcición de Cholesky   

Anteriormente indicamos que el uso de matrices triangulares y la descomposición de Cholesky tienen un papel importante en economía, en R, esta descomposición se obtiene fácilmente M=CholM’ CholM:

MChol<-chol(MVar_Cov)

MChol

##          [,1]      [,2]       [,3]

## [1,] 3.570714 4.0244050 -4.2008403

## [2,] 0.000000 0.7029685  0.5773834

## [3,] 0.000000 0.0000000  0.1398913

t(MChol) %*% MChol

##        [,1]   [,2]  [,3]

## [1,]  12.75  14.37 -15.0

## [2,]  14.37  16.69 -16.5

## [3,] -15.00 -16.50  18.0

Aplicaciones 3. Generar variables aleatorias con autocorrelación

Adicional a la descomposición de la matriz de varianza covarianza de los residuos de los modelos var, para obtener residuos sin autocorrelación, la descomposición de Cholesky se utiliza para genera variables aleatorias con cierta autocorrelación.

# Simula matriz de correlación

C <- matrix(c(1,0.85,0.85,1),2,2)

L <- chol(C)

L

##      [,1]      [,2]

## [1,]    1 0.8500000

## [2,]    0 0.5267827

tau <- diag(c(1,4))

tau

##      [,1] [,2]

## [1,]    1    0

## [2,]    0    4

Lambda <- tau %*% t(L)

Lambda

##      [,1]     [,2]

## [1,]  1.0 0.000000

## [2,]  3.4 2.107131

Z <- cbind(rnorm(5000),rnorm(5000))

cor(Z)

##             [,1]        [,2]

## [1,] 1.000000000 0.009605962

## [2,] 0.009605962 1.000000000

X <- Lambda %*% t(Z)

cor(t(X))

##           [,1]      [,2]

## [1,] 1.0000000 0.8533571

## [2,] 0.8533571 1.0000000

Matriz indempotente

Se dice que A es una matriz idempotente si y sólo si, A*A = A.

Aplicaciones 4. Desviaciones de la media

Este tipo de matriz se utiliza para transformar datos en términos de desviaciones de la media. Es decir, se obtiene una transformación del vector original de datos, de forma tal que cada uno de los nuevos valores obtenidos corresponden al número de unidades alrededor de la media aritmética.

x-mean(x)

## [1] -2  3  1  0 -2

x-I*mean(x)

## [1] -2  3  1  0 -2

x-I*a*t(I)%*%x

## [1] -2  3  1  0 -2

x-a*I*t(I)%*%x

## [1] -2  3  1  0 -2

I*x-a*I*t(I)%*%x

## [1] -2  3  1  0 -2

(I-a*I*t(I))

##      [,1] [,2] [,3] [,4] [,5]
## [1,]  0.8  0.8  0.8  0.8  0.8

M0<-diag(length(I))
M0[M0==0]<- -a    # Propiedad indicada por Greene, p.12
M0[M0==1]<- (I-a*I*t(I))
 M0

##      [,1] [,2] [,3] [,4] [,5]
## [1,]  0.8 -0.2 -0.2 -0.2 -0.2
## [2,] -0.2  0.8 -0.2 -0.2 -0.2
## [3,] -0.2 -0.2  0.8 -0.2 -0.2
## [4,] -0.2 -0.2 -0.2  0.8 -0.2
## [5,] -0.2 -0.2 -0.2 -0.2  0.8

Note que una vez obtenida la matriz idempotente (M0), podemos obtener el vector de desviaciones alrededor de la media, pre-multiplicando el vector original de datos por la matriz obtenida.

round(M0%*%x,0)

##      [,1]
## [1,]   -2
## [2,]    3
## [3,]    1
## [4,]    0
## [5,]   -2

sum(I*x-a*I*t(I)%*%x)

## [1] 0

round(sum(M0%*%x),0)

## [1] 0

x

## [1] 1 6 4 3 1

Suma de desviaciones al cuadrado

La suma de desviaciones de la media al cuadrado:

sum((x-mean(x))^2)

## [1] 18

t(x)%*%M0%*%x

##      [,1]
## [1,]   18

Adicionalmente, se puede obtener la suma del producto de las desviaciones de dos vectores respecto a su media.

y<-c(8,12,16,21,14)

sum((x-mean(x))*(y-mean(y)))

## [1] 8

t(M0%*%x)%*%(M0%*%y)

##      [,1]
## [1,]    8

Ahora, dados dos vectores, se puede construir una matriz de sumas de cuadrados y productos cruzados de desviaciones respecto la media, lo que comúnmente es conocido como matriz de varianzas covarianzas si dividimos entre n o n-1.

mv<-matrix(NA,2,2)
mv

##      [,1] [,2]
## [1,]   NA   NA
## [2,]   NA   NA

mv[1,1]<-t(x)%*%M0%*%x
mv[1,2]<-t(x)%*%M0%*%y
mv[2,1]<-t(y)%*%M0%*%x
mv[2,2]<-t(y)%*%M0%*%y
mv

##      [,1] [,2]
## [1,]   18  8.0
## [2,]    8 92.8

# Matriz Varianza Covarianzas
mv/(length(x)-1)

##      [,1] [,2]
## [1,]  4.5  2.0
## [2,]  2.0 23.2

var(cbind(x,y))

##     x    y
## x 4.5  2.0
## y 2.0 23.2

Determinante de una matriz

La función det(), permite obtener el determinante de una matriz cuadrada.

varcov <- mv/(length(x)-1)
det(varcov)

## [1] 100.4

varcov[1,1]*varcov[2,2]-varcov[2,1]*varcov[2,1]

## [1] 100.4

# Ejemplo 1, Chiang y Wainwright, 2006, p.89.
m4<-cbind(c(10,8), c(4,5))
m4

##      [,1] [,2]
## [1,]   10    4
## [2,]    8    5

det(m4)

## [1] 18

En el caso de una matriz diagonal, el determinante se puede obtener a partir del producto de su diagonal principal (Ecuación 2-47, Greene):

matriz6<-diag(c(1,6,7,9))
matriz6

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    6    0    0
## [3,]    0    0    7    0
## [4,]    0    0    0    9

det(matriz6)

## [1] 378

prod(diag(matriz6))

## [1] 378

Aplicaciones 5. Independencia lineal

Una propiedad interesante del determinante es que cuando las columnas de una matriz son linealmente independientes, su determinante es cero (ver Chiang y Wainwright, 2006).

Recordar que la dependencia lineal entre un conjunto de vectores, se determina cuando existe un conjunto de escalares, tal que, la suma de los productos de dichos vectores es igual a cero. Es decir, dos vectores serán linealmente independientes, si el único conjunto vector de escalar que satisface sum(k*v)=0, es el vector nulo de k (en el ejercicio siguiente k=[3,2]).

# Ejemplo 2, capitulo 4 Chiang y Wainwright, 2006, p.62
v1<-c(2,7) ; v2<-c(1,8) ; v3<-c(4,5)
v1; v2; v3

## [1] 2 7

## [1] 1 8

## [1] 4 5

3*t(v1)-2*t(v2)-v3==0

##      [,1] [,2]
## [1,] TRUE TRUE

Aplicaciones 6.  Optimización de dos variables sin restricciones

Otra aplicación de los determinantes, se obtiene a partir del análisis de su signo. Esto así, para determinar máximos o mínimos en los problemas de optimización, a partir de la prueba de los determinantes para la definición de signo.

En el caso de los problemas de optimización donde tenemos dos o más variables en la función objetivo, luego de que verificar las condiciones de primer orden, las condiciones suficientes de segundo orden se obtienen a partir de los determinantes de la matriz hessiana o matriz de segundas derivadas parciales (simétrica por el teorema de Young) (Ver capítulo 11.4 de Chiang y Wainwright, 2006, p. 312):

Donde los determinantes de las matrices H1,H2,…Hn, se conocen como los menores principales. Dada una matriz hessiana (ver cómo obtener una matrizhessiana en R) tomada de el ejemplo 1 presentado por Chiang y Wainwright (2006), en el libro “Métodos Fundamentales de Economía Matemática”, p.315, dada la siguiente matriz hessiana, podemos obtener los menores principales de la matriz de partir del determinante de las matrices cuadradas de diversos orden que podemos obtener de la matriz hessiana original:

H<-matrix(c(4,1,1,1,8,0,1,0,2), 3, 3, byrow=T)
H

##      [,1] [,2] [,3]
## [1,]    4    1    1
## [2,]    1    8    0
## [3,]    1    0    2

# Obtener los menores principales

H1<-H[1,1]
H2<-det(H[c(1:2),c(1:2)])
H3<-det(H[c(1:3),c(1:3)])

H1; H2; H3

## [1] 4

## [1] 31

## [1] 54

Rango de una matriz

Respecto al uso de los determinantes, se distinguen casos de independencia donde el determinante no se hace cero, por lo que, se utiliza el test de rango para determinar independencia. Esto distingue las matrices de rango completo, con las matrices de rango no completo, dado que el determinante es no nulo, solo si tienen rango completo.

Adicionalmente, note que, en el caso de una matriz, se necesita revisar si cada columna o fila, puede ser expresada como una combinación lineal del resto de columnas. Esta tarea, también suele resumirse a partir de rango de una matriz, que se define como el número máximo de columnas linealmente independientes en dicha matriz.

#Ejercicio 5.1. 5. P.88
matriz7<-matrix(c(7,6,3,3,
                  0,1,2,1,
                  8,0,0,8), 3,4,byrow=T)

matriz7

##      [,1] [,2] [,3] [,4]
## [1,]    7    6    3    3
## [2,]    0    1    2    1
## [3,]    8    0    0    8

qr(matriz7)$rank

## [1] 3

Inversa de una matriz

Dada una matriz cuadrada A, su inversa tiene la propiedad de:

inv(A)A=Ainv(A)=I

Se obtiene como:

inv(A) = 1/|A| adj(A)

En R, la función solve, permite obtener esta inversa, siempre que exista:

solve(varcov)

##             [,1]        [,2]
## [1,]  0.23107570 -0.01992032
## [2,] -0.01992032  0.04482072

Aplicaciones 7.  Resolución de sistema de ecuaciones

Esta matriz inversa, es sumamente útil para resolver sistemas de ecuaciones lineales, imagine el siguiente sistema de ecuaciones, tomado de (Chiang y Wainwright, 2006, p.49):

5a+6b+0.5c=8

12a+7b-1c=4

6a+1b+2c=2

Podemos verificar que este se puede escritor a partir de una matriz de coeficientes (=coef), un vector de variables (=x), y otro vector de constantes (=result): Coef*x=r, por lo que, podemos conocer el vector de incógnitas utilizando la matriz inversa de los coeficientes del sistema: x=coef^(-1)*r.

coef<-matrix(c(5,6,0.5,12,7,-1,6,1,2), 3, 3, byrow=T)
coef

##      [,1] [,2] [,3]
## [1,]    5    6  0.5
## [2,]   12    7 -1.0
## [3,]    6    1  2.0

result<-c(8,4,2)
result

## [1] 8 4 2

#Coef*x=r
# tenemos que: x=coef^-1 r

x<-solve(coef)%*%result

# Nombres a la matriz
rownames(x)<-c("a","b","c")
x

##         [,1]
## a -0.4583333
## b  1.5833333
## c  1.5833333

# Para verificar los resultados anteirores:
coef%*%x

##      [,1]
## [1,]    8
## [2,]    4
## [3,]    2

Singularidad de una matriz

Cuando una matriz tiene inversa, se dice es no singular (Chiang y Wainwright, 2006, p. 86). Por tanto, una matriz debe ser cuadrada, como una condición necesaria de la singularidad. Aunque otra condición suficiente es que los renglones sean linealmente independientes.

H

##      [,1] [,2] [,3]
## [1,]    4    1    1
## [2,]    1    8    0
## [3,]    1    0    2

# Verificar cuadratura
dim(H)

## [1] 3 3

# Rango
qr(H)$rank

## [1] 3

Recordar que cuando una propiedad del determinante, es que cuando una columna es múltiplo de otra, este toma el valor de cero (Chiang y Wainwright, 2006, p. 95). Ahora, la matriz de coeficientes del sistema de ecuaciones, utilizada en el problema de sistema de ecuaciones anterior, debe ser no singular para poder tener solución única (coef*x=d).

Entonces, Coef es no singular y existe una solución única para:

# Determinte
det(H)

## [1] 54

Autovalores y autovectores

El autovalor resulta de desarrollar el determinante de la ecuación característica de una matriz cuadrada determinada |M-λ I|=0. En el caso de los autovalores, la matriz de varianza-covarianza son cuadrada y simétrica, además de positivas semi-definidas. (sus valores propios no son negativos).

En R, se pueden obtener estos valores a partir de la función eigen, que regresa una lista de objetos entre los que se encuentran los autovalores y autovectores, llamados values y vectors, respectivamente.

mv_cov <-mv/(length(x)-1)
mv_cov

##      [,1] [,2]
## [1,]    9  4.0
## [2,]    4 46.4

AutoM9 <- eigen(mv_cov)
AutoM9$values   # AutoM9[[ 1 ]]

## [1] 46.823023  8.576977

AutoM9$vectors  # AutoM9[[ 2 ]]

##           [,1]       [,2]
## [1,] 0.1051692 -0.9944543
## [2,] 0.9944543  0.1051692

Algunas propiedades interesantes de los valores propios son:

•             La suma de los autovalores es igual a la traza de la matriz.

sum(AutoM9$values)

## [1] 55.4

#sum(diag(AutoM9))

•             Lo anterior se aplica al cuadrado:

sum(AutoM9$values^2)

## [1] 2265.96

#sum(diag(AutoM9)^2)

•             El producto de los autovalores, es igual al determinante de la matriz.

det(mv_cov)

## [1] 401.6

prod(AutoM9$values)

## [1] 401.6

•             El Rango es igual al número de autovalores distintos de cero.

qr(mv_cov)$rank

## [1] 2

sum(AutoM9$values != 0)

## [1] 2

Aplicaciones 8. Análisis factorial y componentes principales

Insertando una base de datos hipotética:

n<-50
edad<-runif(n, min = 0, max = 1)
ingresos<-rlnorm(n, meanlog = 10, sdlog = 2)
sexo<-c(rep("Hombre", 25), rep("Mujer", 25))
random_error<-rnorm(n,0,1)
peso<-100+0.8*edad+random_error

data<-data.frame(edad, ingresos, peso)

#Se normalizan las series, para evitar el efecto escala.

data1<-scale(data)

head(data)

##        edad   ingresos      peso
## 1 0.6826406   8018.602  99.58492
## 2 0.5031591  28284.382 100.17860
## 3 0.6639757  44836.382 100.60202
## 4 0.4613249 506105.901 100.03208
## 5 0.5327380   6093.081 100.37986
## 6 0.2899249 111988.762 101.27128

head(data1)

##               edad   ingresos        peso
## [1,]  0.5507769212 -0.4104752 -1.10678996
## [2,] -0.1088741632 -0.3542879 -0.45405375
## [3,]  0.4821776527 -0.3083971  0.01148658
## [4,] -0.2626276643  0.9704823 -0.61515173
## [5,] -0.0001624125 -0.4158137 -0.23277964
## [6,] -0.8925768479 -0.1222158  0.74731739

Posteriormente, se deben obtener los autovalores y autovectores, a partir de la matriz de varianza covarianza de los datos:

MVarVoc<-cov(data1)
MVarVoc

##                edad    ingresos       peso
## edad     1.00000000  0.05681096  0.1964776
## ingresos 0.05681096  1.00000000 -0.1650269
## peso     0.19647760 -0.16502693  1.0000000

Eigenvalues <- eigen(MVarVoc)$values
Eigenvalues

## [1] 1.2303762 1.0558672 0.7137566

Eigenvectors <- eigen(MVarVoc)$vectors
Eigenvectors

##            [,1]        [,2]       [,3]
## [1,]  0.5353602  0.63552274  0.5563275
## [2,] -0.4006233  0.77092776 -0.4951478
## [3,]  0.7435660 -0.04220469 -0.6673293

Note, que estos representan una descomposición de la matriz original:

d <- diag(Eigenvalues)
e <- Eigenvectors

f  <- e %*% d %*% t(e)
f

##            [,1]        [,2]       [,3]
## [1,] 1.00000000  0.05681096  0.1964776
## [2,] 0.05681096  1.00000000 -0.1650269
## [3,] 0.19647760 -0.16502693  1.0000000

Una vez obtenidos los autovalores y autovectores de la matriz de varianza covarianza, o matriz de correlaciones, los componentes principales se obtienen a partir de una multiplicación de matrices.

PC <- as.matrix(data1) %*% Eigenvectors

En tal sentido, los Eigenvectors corresponden a ponderadores de cada variable. Según el valor asumido por dicho ponderador, es posible asignarle interpretaciones a los componentes, por ejemplo, en el caso de los tipos de interes, puede resultar que el primer componente pondere mas los tipos a corto plazo.

pc1 <- data1[,1] * Eigenvectors[1,1] +
       data1[,2] * Eigenvectors[2,1] +
       data1[,3] * Eigenvectors[3,1]

all(PC[,1] == pc1)

## [1] TRUE

Estos componentes no se encuentran autocorrelado.

pairs(PC)

Se verifica, que porcentaje de la varianza total de la matriz original de datos, es explicada por cada uno de los autovalores obtenidos.

PorcentExp<-Eigenvalues/sum(Eigenvalues) * 100
names(PorcentExp)<-c("CP1", "CP2", "CP3")
PorcentExp

##      CP1      CP2      CP3
## 41.01254 35.19557 23.79189

Potencia de una matriz

library(expm)

coef%^%2

##      [,1]  [,2] [,3]
## [1,]  100  72.5 -2.5
## [2,]  138 120.0 -3.0
## [3,]   54  45.0  6.0

Aplicaciones 9. Cadenas de Markov finitas

Los procesos de Markov se emplean para medir o estimar movimientos en el tiempo (Chiang y Wainwright, 2006, p.78). En este sentido, las matrices combinan probabilidades de estado, en las denominadas probabilidades de transición (M), combinado con la distribución de las observaciones en un momento en el tiempo (xt), para ir describiendo la distribución de las observaciones en cualquier momento futuro.

xt’ M = xt+1’

# Ejemplo 1. Chiang y Wainwright, 2006, p.80
# Distribución inicial de empleados

xt0<-c(100,100)
xt0

## [1] 100 100

# Matriz de transición

M<-matrix(c(0.7,0.4,0.3,0.6),2,2)
M

##      [,1] [,2]
## [1,]  0.7  0.3
## [2,]  0.4  0.6

# La Distribución en t=1

xt1 <- xt0%*%M
xt1

##      [,1] [,2]
## [1,]  110   90

# La distribución después de dos periodos
xt2 <- xt0%*%(M%^%2)
xt2

##      [,1] [,2]
## [1,]  113   87

# La distribución después de tres periodos
xt3 <- xt0%*%(M%^%3)
xt3

##       [,1] [,2]
## [1,] 113.9 86.1

Dado lo anterior, se puede verificar que mediante un loops for, podemos obtener la distribución a lo largo de k periodos:

k<-5
for (i in 0:k){
 xt<-xt0%*%(M%^%i)
 print(xt)
}

##      [,1] [,2]
## [1,]  100  100
##      [,1] [,2]
## [1,]  110   90
##      [,1] [,2]
## [1,]  113   87
##       [,1] [,2]
## [1,] 113.9 86.1
##        [,1]  [,2]
## [1,] 114.17 85.83
##         [,1]   [,2]
## [1,] 114.251 85.749

Adicionalmente, podemos estar interesado en la probabilidad de llegar al estado j, después de n pasos (pjn). Por lo que, si tenemos una persona desempleada hoy, sabiendo que si una persona está desempleada actualmente:

La probabilidad de que siga desempleado luego de n periodos es:

Así, las probabilidades de que un individuo que hoy este desempleado, siga desempleado luego de tres periodos:

xt0<-c(0,1)
xt0%*%(M%^%3)

##       [,1]  [,2]
## [1,] 0.556 0.444

Referencias

Chiang, Alpha and Wainwright, Kevin (2006). Métodos fundamentales en economía matemática. 4ª ed. México D.F. México.

Friendly, Michael (2018). Eigenvalues and Eigenvectors: Properties. Recuperado de: https://cran.r-project.org/web/packages/matlib/vignettes/eigen-ex1.html. Consultado en 4/15/2019

Github (2018). Generating correlated random variables with the Cholesky factorization. Consultado el 16.4.2018. Recuperado de: https://mattelisi.github.io/post/generating-correlated-random-variables-with-the-cholesky-factorization/

Greene, William H. (). Análisis econométrico.

Marín, J. y Rubio, G. (2011). Economía financiera. Barcelona: Antoni Bosch.

Sydsaeter, Knut and Hammond, Peter (1996). Matemáticas para el Análisis Económico. Pretince Hall. Madrid.

Tang, Dave. Step by step Principal Component Analysis using R. Recuperado de: https://davetang.org/muse/2012/02/01/step-by-step-principal-components-analysis-using-r/. Consultado en 4/15/2019

Wooldridge, Jeffrey (2010) Introducción a la econometría un enfoque moderno. Cengage Learning. México D.F.

(2009). Principal Component Analysis using R. Recuperado de: http://statmath.wu.ac.at/~hornik/QFS1/principal_component-vignette.pdf. Consultado en 4/15/2019