MATLAB ( Matrix Laboratory: "laboratorio de matrices") es uno de los software matemáticos más utilizados alrededor del mundo y cuenta con su propio lenguaje de programación: el lenguaje m. Esta herramienta es útil para el análisis matemático de variables, la optimización y la simulación. Por ello, este proyecto de estudio busca, además de introducir al usuario en la programación de MATLAB, enseñar el uso de este software para la optimización en economía.

 

El proyecto se divide en 4 partes: la primera de ellas es un pequeño resumen de cómo programar en MATLAB y cómo hacer uso de sus herramientas matemáticas y graficas; la segunda otorga fundamentos en optimización estatica en economía aplicada en el software; la tercer enseña dinámica y como hacer uso de ella para la optimización; por último, la cuarta muestra el uso de MATLAB para resolver problemas de programación dinámica.

Si se desea conocer más sobre las herramientas que ofrece MATLAB, su documentación y sus diferentes aplicaciones, se puede consultar la pagina oficial del software: https://la.mathworks.com/products/matlab.html

 

  

Parte I

 

Programación en MATLAB

 

El lenguaje m de MATLAB es de alto nivel y se encuentra orientado a objetos. Esto quiere decir que se pueden utilizar todas las herramientas tradicionales de la programación para construir diferentes estructuras funcionales. Además, existen funciones ya construidas que facilitan el uso de MATLAB, reduciendo una gran serie de problemas al uso de comandos. Así, el aprendizaje de la sintaxis del lenguaje m es relativamente sencillo.

 

A continuación, se indican algunos fundamentos básicos de la logica de programación en MATLAB y sus comandos más usados para optimizar. 

 

Operadores:

 

Exprensión en MATLAB	Operación
+  y  -	Suma y resta respectivamente
*  y  /	Multiplicación y división
<  y  >	"Menor que" y "mayor que"
=<  y  >=	''Menor o igual que'' y ''mayor o igual que''
==  y  ~=	"Igual a" y "diferente a"
&	Expresión lógica "y" (conjunción)
|	Expresión lógica "o" (disyunción)
~	Negación
'	Transposición
^	Elevar a una potencia

 

Funciones Matemáticas:

 

Comando	Función
Sin / Sinh	Seno / Seno hiperbólico
Cos / Cosh	Coseno / Coseno hiperbólico
Tan / Tanh	Tangente / Tangente hiperbólica
Sec / Sech	Secante / Secante hiperbólica
Csc / Csch	Cosecante / Cosecante hiperbólica
Cot / Coth	Cotangente / Cotangente hiperbólica
Log / Exp	Logaritmo natural / Exponencial
Round /fix	Redondear al entero más cercano / Redondea a cero
Mod / rem	Operador modulo / Residuo de la división
Sqrt	Raiz cuadrada
Floor / Ceil	Redondea hacia abajo / Redondea hacia arriba
diff(x)	Primera diferenciación de x

 

 Además de estas, existen multiples comandos aritmeticos que permiten realizar todo tipo de operaciones matematicas básicas. Estos se pueden consultar por medio del comando Help elfun.

 

Variables

 

La logíca de Matlab funciona por medio de vectores y matrices, por lo que cada variable que se cree será en realidad un vector o arreglo. Para crear una variable se usa el operador de asignación "=", que funciona de derecha a izquierda, por lo que x=7 tiene sentido, pero 7=x no.  Así, Matlab mostrará un vector de tamaño 1x1 con una variable real (double) llamada x que toma el valor de 7.

 

Vectores

 

Existen dos formas para crear vectores en Matlab. La primera de ellas consiste en generar un arreglo con una secuencia de numeros:

 

>> y=1:8

 

Este comando generara un vector y de tamaño 1x8 con los numeros 1,2,3,4,5,6,7 y 8.  Si se desea que la diferencia entre un número y otro sea diferente a uno, es necesario introducir otro argumento intermedio:

 

>> y=1:0.5:8

 

Así, se genera un vector de 1x15 que contiene los números 1, 1.5, 2, 2.5 ,3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5 y 8.

 

La segunda forma para crear un vector consiste en introducir los elementos del vector de la siguiente forma:

 

>> y=[5 -2 39]

 

Este comando crea un vector de tamaño 1x3 con los numeros 5, -2 y 39.

 

Note que hasta ahora solo hemos creado vectores fila. Para generar un vector columna, es necesario crear la transposición de un vector fila ( lo cual se logra con el operador ')  o usar un punto y coma ( " ; ") entre cada elemento del vector.

 

>> z=y' 

>> z= [ 5 ; -2 ; 39]

 

z es un vector columna de tamaño 3x1.

 

Operaciones con vectores 

 

Los vectores en Matlab permiten todas las operaciones tradicionales del algebra lineal. Además, es posible usar una multiplicación elemento a elemento. A continuación, se resumen las operaciones básicas usando los vectores a=[1 2 3] y b=[3 4 5].

 

>> c= a+b = [4 6 8 ]

>> d= a-b = [-2 -2 -2]

 

Para sumar y restar es necesario que los vectores a y b tengan el mismo tamaño. Por otro lado, cuando se desea realizar una multiplicación, se deben respetar las normas de la multiplicación matricial. Así, a*b no es posible realizarla, mientras a*b' sí.

 

>> e= a*b' = 26

 

Si se desean realizar multiplicaciones entre vectores elemento a elemento, basta con escribir un punto ( " . ") despues de las variables, así:  

 

>> f= a.*b = [3 8 15]

 

Esto es útil cuando se desea encontrar la potencia de un vector elemento a elemento, pues basta con usar a.^2. De la misma forma, la división elemento a elemento se genera con el punto.

 

>> g= a./b= [0.33 0.50 0.60]

 

Para indexar un vector, es decir, acceder a una posición o a una elemento en una posición, basta con escribir el nombre del vector, usar los corchetes "( )" y dentro de ellos escribir el indice o posición:

 

>> b(2)

 

El comando imprimiría 4, que es el segundo elemento de b. Esto sucede porque Matlab cuenta los indices empezando en 1, a diferencia de la mayoria de lenguajes de programación donde inician en 0.

 

Para modificar un elemento de un vector, basta con asignarle un valor a la indexación de ese vector en alguna posición:

 

>> b(2)=-7

 

Así, b ahora es b=[3 -7  5].

 

Por último, los vectores admiten operaciones con escalares:

 

>> a*2

>> a+2

>> a-2

>> a/2

 

 

Matrices

 

Para crear una matriz, se debe generar una sucesión de vectores separados por un punto y coma ( " ; "). Por ejemplo, para crear una matriz "A" de tamaño 3x3 usariamos:

 

>> A= [-2 3 1 ; 3 7 21 ; -65 102 21]

 

 La matriz A es transponible por medio del apostrofe " ' ".

 

>> B=A'

 

Además, también pueden sumarse, restarse o multiplicarse matrices:

 

>> A+B

>> A-B

>> A*B

 

Para acceder a elementos de la matriz, debe usarse el nombre de la matriz y entre parentesis señalarse el numero de fila ( " i ") en la que se encuentra el elemento requerido y su columna ( " j " ).  Así, la estructura es A(i, j). Para acceder al elemento de la segunda fila y primera columna, por ejemplo, se usa:

 

>> A(2, 1)

 

Y para modificarlo, basta con igualar el llamado a un valor:

 

>> A(2, 1) = 5

 

Para acceder a toda una fila, se puede usar el operador dos puntos ( " : ") de la siguiente forma:

 

>> A(1, : )

 

De forma similar, si se quiere retornar una columna:

 

>> A( :, 1)

 

 

Comandos útiles en vectores y matrices

 

Comando	Salida
Length (x)	Largo de un vector x
Det (A)	Determinante de la matriz A
Size (A)	Tamaño de una matriz A
Inv (A)	Inversa de la matriz A
Zeros (n)	Vector de ceros de tamaño n
Zeros (m,n)	Matriz de ceros de tamaño mxn
Eye (n)	Matriz identidad de tamaño nxn
Rand	Un vector de números seudoaleatorios entre 1 y 0
Rand (m,n)	Una matriz de números seudoaleatorios entre 1 y 0 de tamaño mxn

 

Graficar

 

 Las gráficas permiten observar el comportamiento de las funciones, lo que ayuda a intuir algunas características de la función, tales como sus puntos criticos, su forma o su convexidad. Para crear una gráfica en 2d, es necesario hacer uso del comando plot().

 

 >> x= rand(20,1)

 >> plot(x)

 

En ocasiones es necesario crear más de una gráfica en una misma figura, para ello se debe usar el comando Hold on y luego hacer el plot, por ejemplo:

 

>> hold on

>> y= rand(20,1)

>> plot(y)

 

Ahora se mostrarán ambas gráficas en una sola figura.

 

Podemos agregar a la figura varios tipos de especificaciones. Entre ellos, de linea, titulos y colores. Por ejemplo, si se desea que el plot de x sea una línea descontinua azul hacemos:

 

>>plot(x,'b--')

 

Donde b es por "blue" ( azul en inglés) y "--" representan una linea discontinua. 

Para agregar titúlos, títulos en los ejes y leyendas se usa:

 

>> title('escriba aquí el título')

>> xlabel ('escriba aquí el título del eje x')

>> ylabel('escriba aquí el título del eje y')

 

Funciones

 

Para definir una función se utilizan los editores de texto o "scripts".  Un script permite realizar un programa que puede guardarse y ejecutarse en su totalidad y no línea por línea como con la consola.  Crear una función con un script es sencillo, debe abrirse un nuevo script y escribir:

 

>> Function f=ejemplo01(x, y)

>> f=x.^(1/2) + y.^(1/2)

>> end

 

El comando "function" crea una función "f" que puede ser llamada a través del nombre "ejemplo01" y toma como argumentos a "x" y a "y". Es importante que al guardar el script, el nombre de la función y del archivo sean los mismos.

 

Para imprimir valores de una función, es necesario llamarla en otro script con su nombre. En este caso: "ejemplo01". 

 

 a=-30:1:30;

 b=a;

  q=ejemplo01(a,b);

 

 Así, el valor de q correspondrá a la evaluación de a y b en la función que definimos en el script "ejemplo01".

 

Si queremos realizar una gráfica de q, debido a que estaría en 3D, es necesario utilizar "meshgrid" que genera una cuadrícula con los datos de dos arreglos x y y. La cuadrícula o rejilla tendrá lenght(x) divisiones en el eje x y lenght(y) divisiones en el eje y. Luego, a cada punto de la rejilla se le asigna un valor de la funcion "ejemplo01". Por último, se hace uso del comando "mesh" que genera el gráfico a partir de las cuadriculas. 

 

 x=-30:1:30;

 y=x;

 [X,Y]= meshgrid(x,y);

F=ejemplo01(X,Y);

 mesh(X,Y,F)

 

Se mostrará en una figura la gráfica de la función como se definió previamente.  Si se quieren simular como cambia la función cuando cambian los valores de los parametros, basta con cambiar los parametros a un valor cualquiera "a" y utilizar un ciclo for en el script principal. Así en "ejemplo01" tendriamos:

 

 

 function f=ejemplo01(x,y,a)

 f=x.^a + y.^(1-a)

 

 

En el segundo script:

 

 x=0:1:20;

 y=x;

 a=0:0.1:10

 [X,Y]=meshgrid(x,y)

 For i=1:lenght(a)

 F=ejemplo01(x,y,a(i));

 mesh(X,Y,F);

 pause(0.5);

 end

 

 Otra forma de definir funciones es por medio del arroba "@". Este comando permite definir una función en un script y utilizarla en el mismo. Para implementarla basta con poner el simbolo "@" y, después, entre parentesis, los argumentos que toma la función. Por último, se agrega la forma funcional que se desee. A continuación  se muestra cómo realizar este proceso mediante un ejemplo:

 

Abrimos un script y escribimos:

 

fun= @(x) 100*( x(2) - x(1)^2) ^2 + (1-x(1))^2;

 

Luego podemos llamar a la función utilizando el nombre "fun".

 

Parte II

 

Métodos de optimización estática

 

La optimización estática es un conjunto de métodos análiticos y númericos que buscan encontrar el mejor valor para una o varias variables, de tal forma que estos generen un valor óptimo en un problema. El valor óptimo será aquel que resuelva el problema de optimización, que  en economía, generalmente, es un problema de decisión. 

En Matlab existen varias rutinas que, por medio de métodos númericos, optimizan una función. Entre ellos encontramos: 

 

 1) fminbnd: Rutina que encuentra el mínimo de una función en una variable por medio de  un proceso de optimización no restringida. fminbnd toma principalmente tres argumentos: fx, X₁ y X₂. El primero se refiere a la función que se quiere optimizar, mientras X₁  y  X₂ son constantes que forman un intervalo sobre el cual se optimiza. Así, fminbnd devuelve un argumento minizador  X tal que X₁<  X < X₂.

2) fminsearch: Rutina que calcula el mínimo de una función de varias variables sin restricciones.  Su forma funcional es similar a la de fminbnd. 

 3) fmincon: Rutina que haya los argumentos minimizadores de una función en varias variables sujeta a unas restricciones de la forma AX =< b. Donde A es una matriz de coeficientes, x un vector de variables y b un vector de constantes. 

 

A continuación se muestran tres ejemplos que ejecutan las rutinas expuestas anteriormente:

 

1) En un nuevo script creamos una función llamada bnd:

function f= bnd(x)

f=x.^3 + 5*x.^2 +2;

end

 

En un segundo script creamos los datos:

x=[-5:0.01:-1];

 

Se aplica la función:

 

f=bnd(x);

 

Se haya el mínimo:

 

[xmin, fmin,ef,pt]= fminbnd('ej6',-1.5,-0.5)

 

Como se puede observar, las salidas que arroja la rutina fminbnd son cuatro. La primera de ellas contiene el argumento minimizador. La segunda, por su parte, se refiere al valor minímo de la función f. La tercera (ef)  es una constante tal que si ef > 0, el algoritmo de solución convergió a la solución. Esto quiere decir que se hayó la solución exacta del problema la iterar el método númerico sobre el problema. Si ef<0, la solución no convergió y si ef=0, no se pudo encontrar una solución tal que el método convergiera a ella.  La cuarta salida (pt) otorga información extra sobre el método usado y la solución. 

 

 En cuanto a los argumentos de fminbnd, estos son la función, seguida de los valores X1 y X2 a los que nos referimos anteriormente.

2) Para esta rutina utilizaremos la segunda forma de definir funciones y que introducimos anteriormente. Esta es en la que usamos el arroba " @ ".

 

En un nuevo script escribimos: 

 

func= @(x) 100*( x(2) - x(1)^2) ^2 + (1-x(1))^2;

 

 A una variable X0 le asignamos los valores X1 y X2:

 

x0 = [-1.2, 1 ];

 

Finalmente hallamos el óptimo. Para esto llamamos a la función func como un argumento de fminsearch usando las comillas simples, luego colocamos el rango sobre el que se minimiza, dado por x0. 

 

[x,f]= fminsearch('func',x0);

 

Las salidas son x y f. La primera es el argumento minimizador y la segunda el valor de la función.

 

3)  Note que para utilizar una rutina de optimización basta con tres cosas: la función, los valores iniciales y la rutina con sus respectivas entradas y salidas. Este proceso no varía cuando se usa fmincon, solo que en lugar de crear los datos de X, tambien deben crearse los de A y b:

 

func= @(x) 100*( x(2) - x(1)^2) ^2 + (1-x(1))^2;

 

x=[-1,2];

A=[1,2];

b=1;

 

xmin=fmincon('func',x,A,b);

 

La salida xmin es el argumento minimizador. fmincon puede arrojar varias interesantes salidas, tales como: 

 

Salida	Descripción
xmin	Argumento minimizador
fval	Valor mínimo de la función
exitflag	Describe la condición de salida de fmincon
output	Información sobre el proceso de optimización
lambda	Multiplicadores de lagrange o, como se les conoce en economía, precios sombra
grad	Gradiente de la función en la solución
hessian	Matriz hessiana

 

 

Parte III

 

Ecuaciones Diferenciales Ordinarias (ODES)

 

 Matlab posee multiples rutinas para resolver ODES, entre ellas se encuentran ode45 y ode23tb. Ode45 es un método de solución de ecuaciones diferenciales ordinarias y su sintaxis es:

 

[t,x]=ode45(@función,tinit,xinit)

 

Donde @función es la función que contiene la ecuación diferencial a resolver, tinit es el intervalo de tiempo sobre el que se resuelve la ecuación y xinit es un vector que actua como condiciones iniciales para el problema.  Las salidas t y x corresponden, respectivamente, a los valores del tiempo y la solución de la función. 

 

Es importante recordar que la solución de una ecuación diferencial es una función del tiempo, es decir, una trayectoria tal que para cada valor de t, existe un valor de la función. 

 

 Con el fin de configurar más salidas, es posible usar options. Options es configurada por medio del comando odeset, que genera salidas tales como la jacobiana, gráficos y estádisticas.  Un ejemplo de esto es:

 

options=odeset('outputfcn','odeplot')

[t,x]=ode45(@función,tinit,xinit, options)

 

Options hará que se arroje un gráfico de la solución al problema.  Las opciones son muy amplias y la explicación sobre ellas en la página oficial de Matlab es  detallada, por lo que no resulta útil resumirlas aquí.

 

A continuación se muestra un ejemplo del uso de Ode45 para resolver un problema:

 

En un script escribimos la función:

 

function dxdt=deq(t,x)

dxdt=3-x

end

 

En un segundo script:

 

tinit=[0 20]

xinit=[-1.5];

options=odeset('outputfcn','odeplot')

[t,x]=ode45(@deq,tinit,xinit,options)

 

La solución a los problemas de ecuaciones diferenciales varía con respecto a las condiciones iniciales y a los parametros que posea la función.  Por esto, a continuación se muestra una forma simple de incluir un cambio en los parametros en la solución. Para las condiciones iniciales el proceso es similar, por lo que no se explica. 

 

En un primer script definimos la función:

 

 

Sistemas de Ecuaciones diferenciales Ordinarias

 

Anteriormente vimos el uso de Ode45 para resolver una ecuación diferencial de primer orden. Sin embargo, esta rutina permite dar solución a ecuaciones hasta de orden n. Para esto, se debe transformar la ecuación en un sistema de ecuaciones ordinarias de primer orden y luego implementarlo en Matlab. A continuación se muestra un ejemplo de como hacer esto:

 

Dada la ecuación:

 

 

Se sustituye:

 

 

  

 

 

  

Note que son tres ecuaciones porque la variable y se encontraba derivada tres veces. Para implementar este sistema en Matlabs se crean tres ecuaciones con los terminos del sistema dentro de una función y agrega una más que representa un vector de tres filas por una columna, como se muestra a continuacón:

 

En un nuevo script:

function dydt=sist(t,x)

dydt=zeros(3,1)

dydt=x(1)

dydt=x(2)

dydt=1-4^3 - 2*x(3) + x(2) + 2*x(1)

end

 

En general, el sistema tiene tantas ecuaciones como se encuentre derivada la variable de interes y el tamaño del vector de ceros es también determinado por esta cantidad.  Para darle solución, se sigue el mismo procedimiento que con las ODES. Sin embargo, hay que tener cuidado al defini xinit, pues esta debe tener tantas condiciones iniciales como ecuaciones tiene el sistema. 

 

En un segundo script:

tinit=0:0.5:10

xinit=[1 0 0]

[t,x]=ode45(@sist,tinit,xinit)

 

Diagramas de fase

 

Para realizar un diagrama de fase de un sistema de ecuaciones diferenciales ordinarias, basta con realizar un ciclo que genere y grafique las soluciones del sistema para una serie de condiciones iniciales. Así, se observa si las soluciones convergen, divergen o son del tipo lyapunov, lo que permite hacer el análisis cualitativo del sistema.  Dado esto, basta con modificar el anterior script de la siguiente forma:

 

tinit=0:0.5:10;

xinit=linspace(-1,1,10);

figure;

hold on

 for i=1:10

for j=1:10

for h=1:10

[t,x]= ode45(@sist,tinit,[xinit(i) xinit(j) xinit(h])

 plot(x(:,1),x(:,2),x(:,3))

end

end

end

 

La figura deberá contener un diagrama de fase del sistema de ecuaciones que se planteó anteriormente. 

  

Parte IV

 

 

Introducción a la programación dinámica

 

Con el fín de conocer como se resuelve un problema de programación dinámica en Matlab, se seguirá la metodología planteada por Seguí, M (sin fecha), quien resume como hacer uso de métodos numéricos para la solución de diferentes problemas de optimización. Entre ellos, el planteamiento de la ecuación de Bellman y el uso de Matlab para solucionar el problema del planificador social. A continuación, se encuentra una breve descripción de la ecuación de Bellman, también llamada la ecuación de la programación dinámica. 

 

La ecuación de Bellman y los modelos dinámicos

 

Un modelo dinámico es tal que se encuentra representado por la siguiente estructura de un problema de optimización:

 

 

 Donde f es la función objetivo que depende de una variable de estado x, que representa como cambia el problema a través del tiempo; una variable de control u , que es la variable de decisión, y el tiempo t. Por su parte, beta es una tasa de redescuente, que trae a valor presente los valores de la función en el tiempo k. La sumatoria se usa  para tomar todos los periodos de tiempo en el problema y, como el problema es discreto, es una sumatoria y no una integral. Por último, la primera restricción muestra como cambia el estado y la segunda es una condición inicial para el problema. 

 

Note que beta debe ser menor a uno para que haya convergencia. 

 

La solución a este problema esta dada na función valor. La función valor V(x₀) es el valor maximo del problema anteriormente planteado. Como se puede apreciar, el problema esta expresado en infinitos tiempos, lo que hace que resolver el problema sea muy complicado. Por esto, Bellman introdujo una condición necesaria que simplifica el problema a uno de dos tiempos. 

 

 

Esta ecuación quiere decir que, si existe una trayectoria óptima para la solución del problema y se ha elegido estano importa desde donde se solucione el problema, la trayectoria óptima sigue siendo la misma. En otras palabras, como indica Palomino, P(2014): "Una política óptima tiene la propiedad de que, cualesquiera sean el estado y las decisiones iniciales tomadas, las restantes decisiones deben constituir una política óptima con independencia del estado resultante de la primera decisión". Basandose en esto, Bellman reduce el problema a través de una definición recursiva de la función valor. Ahora basta con optimizar la ecuación de Bellman para el problema que se posea para llegar a una solución, que será una función en el tiempo, la cual indica la decisión a tomar en cada momento t.

  

A continuación se muestra un ejemplo desarrollado por Seguí,M(sin fecha) sobre el uso de la ecuación de Bellman en matlab para resolver el problema del planificador social en una economía. Para conocer todas las condiciones e implicaciones del problema del planificador social, consulte el paper de Seguí, M ( sin fecha) donde se desarrolla todo el problema y se explican las implicaciones económicas de este.

 

 Se crea un archivo llamado vfiterm.m que contendrá el programa para resolver el problema.

 

La ecuación de Bellman a resolver es:

 

 v(k) = max fln(A*k^alpha+(1-delta)*k-k') + beta v(k')g

 con k0 dado

 

Se definen los parametros:

 

 

alpha = 0.35;

beta = 0.95;

delta = 0.06;

A = 10;

 

Luego los parametros extra, como el tamaño del grid, el número de iteraciones, un criterio de tolerancia y el número de periodos a simular.

maxit = 1000;

p = 500;

crit = 1e-2;

T = 100;

Luego se define el capital en estado estacionario kss y las condiciones de inicio, mínimo y máximo:

 

kss = ((A*beta*alpha)/(1-(1-delta)*beta))^(1/(1-alpha));

k0 = (1/2)*kss;;

Kmin = 0.5*kss;

Kmax = 1.2*kss;

K(1)=k0;

K = [Kmin+(0:p-2)*(Kmax-Kmin)/(p-1) Kmax]';

 

Luego se coloca un punto de partida para la función valor y se crea una matriz M:

 

V0 = 0*K;

e = ones(1,p);

M = log(max((A*K.^alpha+(1-delta)*K)*e-e'*K',1e-8));

 

Se itera la función valor como método númerico:

 

 

figure(1)

iconv = 0;

it = 1;

while (iconv==0 & it<maxit)

[V1,G] = max((M+beta*(V0*e)')');

V1 = V1';

G = G';

if norm(V0-V1)<crit

iconv = 1;

end

plot(K,V1)

axis([Kmin Kmax 0 100])

title(['Función de valor en la iteración ',int2str(it),' con error ',num2str(norm(V0-V1))])

xlabel('k')

ylabel('v(k)')

grid

pause(0.01)

V0 = V1;

it = it+1;

end

figure(2)

plot(K,K(G),K,K,':')

title('Regla de decisión óptima')

xlabel('k')

ylabel('kp=g(k)')

grid

text(kss,kss,'o kss')

 

Se obtienen los estados estacionarios de producto, consumo e inversión, luego se define la tasa de ahorro:

 

yss=A*kss^alpha

css=yss+[(1-delta)-1]*kss

iss=yss-css

sss=iss/yss 

 

Se simulan las trayectorias que estos siguen:

 

figure(3)

kt = zeros(T+1,1);

ind = zeros(T,1);

[aux,ind(1)] = min(K<k0);

kt(1) = K(ind(1));

for t=1:T

ind(t+1) = G(ind(t));

kt(t+1) = K(ind(t+1));

end

kss2=kss*ones(T,1)';

z=[1,T,K(1),kss+1];

plot(1:1:T,kt(1:T),'.-',1:1:T,kss2,'r')

title('Trayectorias óptimas para el Capital')

xlabel('t')

ylabel('k t')

grid

axis(z)

legend('Senda', 'SS')

yss2=yss*ones(T-1,1)';

figure(4)

yt = A*kt(1:T).^alpha;

plot(1:1:T,yt,'.-',1:1:T-1,yss2,'r')

title('Trayectorias óptimas para el Producto')

xlabel('t')

ylabel('y t')

grid

legend('Senda', 'SS')

iss2=iss*ones(T-1,1)';

figure(5)

it = kt(2:T+1)-(1-delta)*kt(1:T);

plot(1:1:T,it,'.-',1:1:T-1,iss2,'r')

title('Trayectorias óptimas para la inversión')

xlabel('t')

ylabel('it')

grid

legend('Senda', 'SS')

css2=css*ones(T-1,1)';

figure(6)

ct = yt-it;

plot(1:1:T,ct,'.-',1:1:T-1,css2,'r')

title('Trayectorias óptimas para el consumo')

xlabel('t')

ylabel('c t')

grid

legend('Senda', 'SS')

 

Por último, se guardan los resultados principales en un archivo llamado vfout.mat

save E:/vfoutm V0 G kt ct it yt

  

Al correr el programa (para esto basta con escribir el nombre de este en la consola), se mostrarán las figuras creadas que contendrn las trayectorias óptimas de las variables de interés, es decir, producto, consumo, inversion y el capital. Además se mostrará la regla de decisión óptima y la iteración de la función valor. 

 

  

Referencias Bibliográficas

 

Faruk,N(sin fecha). A brief introduction to using ode45 in Matlab. Universidad de Berkley, California. En línea. Recuperado de: http://www.eng.auburn.edu/~tplacek/courses/3600/ode45berkley.pdf

 

Gonzalez, N. Quintero, L. Maldonado, F. Sanchez, E (2009). Economía matemática en Matlab. Universidad nacional de Colombia. Colombia.

 

Katznelson, Y (2015). Solving systems of first order equations with ode45. En línea. Recuperado de: https://www3.nd.edu/~nancy/Math20750/Demos/3dplots/dim3system.html

  

Seguí,M (sin fecha). Métodos numéricos en Matlab. En línea. Recuperado de:  http://ciep.itam.mx/~msegui/notasmishelle.pdf

 

Lynch, S (2003). Dynamical Systems using Matlab. Birkhauser, Bosen, Berlín.