/*
THIS IS A GAUSS PROGRAM ON HOW TO ITERATE ON THE BELLMAN EQUATION
The simple ECONOMIC model:
V(k) = max (ln (AK^alpha - K')) + Beta V(K')
NOTE THAT:
100 % Capital depreciation (delta == 1) There's no delta here
Written by: O. Mikhail (2000)
omikhail@bus.ucf.edu
omikhail@hotmail.com
*/
new;
@ open the output file @
output file = vkmodel.out reset;
@ format command is used to control the way output is printed @
format /M1 /ROS 10,7 ;
@ Stamp the date and time @
print ;
print "Date : " datestr(0);
print "Time : " timestr(0);
print ;
@ Initialize the coefficients @
@ Calibrate the parameters @
alpha = 0.7 ;
beta = 0.96 ;
A = (10^(1-alpha))/(alpha*beta) ;
@ Initialize the vectors and the matrices @
@ k is for Capital @
k = zeros(5,1) ;
@ kp is for Capital prime @
kp = zeros(5,1) ;
@ v is a table (matrix) to hold the value of the value function @
v = zeros(5,5) ;
@ diff is used to capture the difference between vold and vmax @
diff = zeros(5,1) ;
@ Initialize vold with the steady state value (VKSS in Fortran) @
let vold[5,1] = 39.63350 ;
@ vmax will pick the maximum of v over kp for a given k @
@ then use it to compute the next iteration of v @
vmax = zeros(5,1) ;
poskp = zeros(5,1) ;
let supnorm[5,1] = 0.0001 ;
@ iter is the iteration number V0,V1,V2 ... etc @
iter = 1 ;
do while iter <= 20 ;
	@ the loop i and the loop j will create the table @
	@ i is the index for old capital called k (column) @
	i = 1 ;
	do while i <= 5 ;
	k[i]  = 9.97 + (i / 100)  ;
		@ j is the index for the new capital called capital prime kp @
		  j = 1 ;
		  do while j <= 5 ;
		  kp[j] = 9.97 + (j / 100) ;
		@ Check if Consumption is negative @
		@ If it is negative, just put zero in the table @
                if (A*(k[i]^alpha)-kp[j]) <= 0 ;
                v[i,j] = 0 ;
                break ;
                endif ;
		  @ otherwise @
	  	  @ Compute the value function @
		  v[i,j] = ln(A*(k[i]^alpha)-kp[j]) + (beta * vold[j]) ;
		  j = j + 1 ;
		  endo ;
	@ get the maximum of the row i in the table @
	@ this is the max value for kp given k @
  	vmax[i] = maxc(v[i,.]') ;
	@ get the column position of this maximum @
	@ i.e. to get later the value of kp at which this max occur @
	poskp[i] = maxindc(v[i,.]') ;
	@ continue looping over the rows @
	i = i + 1 ;
	endo ;
@ do some fancy print, not necessary, can be commented out @
print "***********************************************************" ;
print " Iter " iter ;
print ;
print " k and kp " ;
print k~kp ;
print ;
print " The V matrix " ;
print v ;
print ;
print " Vmax " ;
print vmax ;
print ;
print " Position " ;
print poskp ;
print ;

@ compute the difference between vmax (new) and vold == diff @
diff = abs(vmax - vold) ;
@ if diff is small ==> convergence , break out of the iter loop @
@ z will equal 1 if all the elements of diff are le (less or equal) supnorm @
z = dotfle(diff,supnorm) ;
print " z diff supnorm ";
print z~diff~supnorm ;
print ;
	if z == 1 ;
		print " ************************************* " ;
		print " *********** Got you ***************** " ;
		print " ************************************* " ;
		print ;
		print " ITERATION " iter ;
		print ;	
		print " Vold and Vmax " vold~vmax ;
		print ;
		print " ************************************* " ;
		print " Now I am getting out the iter loop " ;
		print " reason : I am really dizzy " ;
		print " ************************************* " ;
		print ;
		break ;
	endif ;

@ Since diff is not small ==> no convergence, then vold == vmax @
print " Vold before and Vmax " vold~vmax ;
print ;
vold = vmax ;
print " Now vold is = vmax " vold~vmax ;
print ;
print ;
@ Repeat the iterations @
iter = iter + 1 ;
endo ;


print "***********************************************************" ;
print ;
print " Iteration number " iter ;
print ;
print " Capital     Capital Prime   Value function ";
print k~(9.97+(poskp/100))~vmax ;
print ;
print "***********************************************************" ;
@ Nice to meet you, please come again @
@ O. Mikhail @
end ;