Merge pull request jstac#2 from tomohitosan/master

Update

edtc_matlab_rev/README.txt

@@ -1,3 +1,6 @@
 The following files are left untouched (i.e. no modification) because they are written in a standard way for Matlab coding.
+- ar1.m
+- ecdf.m
 - fphamilton.m
-- quadmap1.m
+- polyclass0.m
+- quadmap1.m

edtc_matlab_rev/ar1.m

@@ -0,0 +1,12 @@
+% Filename: ar1.m
+% Author: Andy Qi
+% Date: December 2008
+% Corresponds to: Listing 8.1
+
+a = 0.5;
+b = 1;
+X = zeros(1, 101);    % Create an empty array to store path
+X(1) = normrnd(0, 1); % X_0 has dist N(0, 1)
+for t = 1:100
+    X(t + 1) = normrnd(a * X(t) + b, 1);
+end

edtc_matlab_rev/cpdynam_rev1/README.txt

@@ -0,0 +1,3 @@
+Listing 2.2 is divided into two functions, evaluate and differentiate.
+def __init__(self,coef) in the python code is omitted. 
+Instead, "coef" is just dealt with as an input in those two functions.

edtc_matlab_rev/cpdynam_rev1/T.m

@@ -0,0 +1,22 @@
+% Filename: T.m
+% Author: Tomohito Okabe
+% Date: April 2013
+% Corresponds to: T in Listing 6.7
+
+
+function t = T(p, x, W, P, D)
+
+global alpha
+
+% Computes Tp(x), where T is the pricing functional operator.
+% Parameters : p is an instance of lininterp and x is a number.
+
+y = alpha * mean(p(W));
+            if y <= P(x)
+                t = P(x);
+                return;
+            end
+h = @(r) alpha * mean(p(alpha * (x - D(r)) + W));
+t = fix_point(h, P(x), y);
+end
+

edtc_matlab_rev/cpdynam_rev1/cpdynam_rev1.m

@@ -0,0 +1,46 @@
+% Filename: cpdynam_rev1.m
+% Author: Tomohito Okabe
+% Date: April 2013
+% Corresponds to: Listing 6.7
+
+% Note: All of the code is wrapped in a class definition so that it can
+% be put in one file.  See the comments in kurtzbellman.m.  An example
+% of usage is given below.
+
+global alpha a c
+    alpha = 0.8;
+    a = 5.0;
+    c = 2.0;
+
+ W = betarnd(5, 5, 1, 1000) * c + a; % Shock obs.
+ P = @(x) 1.0/x;     % Inverse demand function
+ D = P;
+
+
+% The following test shows how to use the codes. lininterp_rev1.m must 
+% be in the current folder.
+
+% gridsize = 150;
+% grid = linspace(a, 35, gridsize);
+% tol = 0.0005;
+% vals = zeros(1, gridsize);
+% new_vals = zeros(1, gridsize);
+% for i = 1:gridsize
+%     vals(i) = P(grid(i));
+% end
+% while 1
+%     hold on;
+%     plot(grid, vals);
+%     p = @(z)lininterp_rev1(grid, vals,z);
+%     f = @(x) T(p, x, W, P, D);
+%     for i = 1:gridsize
+%         new_vals(i) = f(grid(i));
+%     end
+%     if max(abs(new_vals - vals)) < tol
+%         break
+%     end
+%     vals = new_vals;
+% end
+% hold off;
+
+

edtc_matlab_rev/cpdynam_rev1/fix_point.m

@@ -0,0 +1,14 @@
+% Filename: fix_point.m
+% Author: Tomohito Okabe
+% Date: April 2013
+% Corresponds to: fixpoint in Listing 6.7
+
+function fp = fix_point(h, lower, upper)
+% Computes the fixed point of h on [upper,lower] using fzero, 
+% which finds the  zeros (roots) of a univariate function.
+% Inputs : h is a function and lower and upper are numbers 
+% (floats or integers).
+
+fp = fzero(@(x) x - h(x), [lower, upper]);
+
+end

edtc_matlab_rev/cpdynam_rev1/lininterp_rev1.m

@@ -0,0 +1,17 @@
+% Filename: lininterp_rev1.m
+% Author: Tomohito Okabe
+% Date: April 2013
+% Corresponds to: Listing 6.4
+
+function g = lininterp_rev1(X,Y,z)
+            % Uses MATLAB's interp1 function for interpolation.
+            % The z values are truncated so that they lie inside
+            % the grid points.  The effect is that evaluation of
+            % a point to the left of X(1) returns Y(1), while 
+            % evaluation of a point to the right of X(end) returns
+            % Y(end)
+
+            z = max(z, X(1));
+            z = min(z, X(end));
+            g = interp1(X, Y, z);
+end

edtc_matlab_rev/polyclass0.m

@@ -0,0 +1,7 @@
+
+% Filename: polyclass0.m
+%
+% The MATLAB version of the Python code in polyclass0.py is omitted, since
+% polyclass0.py is not real Python, but rather for illustration purposes
+% only.  Please see polyclass.m
+%