求随机动态规划解法的matlab实例，拜托了

求随机动态规划解法的matlab实例，拜托大家了。

找了一个程序，希望对你有所帮助

1. 
   
2. clear all;
   
3. close all;
   
4. % Define the structural parameters of the model,
   
5. % i.e. policy invariant preference and technology parameters
   
6. % alpha : capital's share of output
   
7. % beta  : time discount factor
   
8. % delta : depreciation rate
   
9. % sigma : risk-aversion parameter, also intertemp. subst. param.
   
10. alpha = .35;
    
11. beta  = .98;
    
12. delta = .025;
    
13. sigma = 2;
    
14. 
    
15. % Number of exogenous states
    
16. gz = 2;
    
17. % Values of z
    
18. z = [4.95
    
19.      5.05];
    
20. % Probabilites for z     
    
21. prh = .5;
    
22. prl = 1-prh;     
    
23. % Expected value of z
    
24. zbar = prh*z(1) + prl*z(2);
    
25. 
    
26. % Find the steady-state level of capital as a function of
    
27. %   the structural parameters
    
28. kstar = ((1/beta - 1 + delta)/(alpha*zbar))^(1/(alpha-1));
    
29. % Define the number of discrete values k can take
    
30. gk = 101;
    
31. k  = linspace(0.90*kstar,1.10*kstar,gk);
    
32. 
    
33. % Compute a (gk x gk x gz) dimensional consumption matrix c
    
34. %   for all the (gk x gk x gz) values of k_t, k_t+1 and z_t.
    
35. for h = 1 : gk
    
36.    for i = 1 : gk
    
37.        for j = 1 : gz
    
38.             c(h,i,j) = z(j)*k(h)^alpha + (1-delta)*k(h) - k(i);
    
39.             if c(h,i,j) < 0
    
40.                 c(h,i,j) = 0;
    
41.             end
    
42.             % h is the counter for the endogenous state variable k_t
    
43.             % i is the counter for the control variable k_t+1
    
44.             % j is the counter for the exogenous state variable z_t
    
45.        end
    
46.    end
    
47. end
    
48. % Compute a (gk x gk x gz) dimensional consumption matrix u
    
49. %   for all the (gk x gk x gz) values of k_t, k_t+1 and z_t.
    
50. for h = 1 : gk
    
51.    for i = 1 : gk
    
52.        for j = 1 : gz
    
53.             if sigma == 1
    
54.                 u(h,i,j) = log(c(h,i,j))
    
55.             else
    
56.                 u(h,i,j) = (c(h,i,j)^(1-sigma) - 1)/(1-sigma);   
    
57.             end
    
58.             % h is the counter for the endogenous state variable k_t
    
59.             % i is the counter for the control variable k_t+1
    
60.             % j is the counter for the exogenous state variable z_t
    
61.        end
    
62.    end
    
63. end
    
64. 
    
65. % Define the initial matrix v as a matrix of zeros (could be anything)
    
66. v = zeros(gk,gz);
    
67. % Set parameters for the loop
    
68. convcrit = 1E-6;  % chosen convergence criterion
    
69. diff = 1;         % arbitrary initial value greater than convcrit
    
70. iter = 0;         % iterations counter
    
71. while diff > convcrit
    
72.     % for each combination of k_t and gamma_t
    
73.     %   find the k_t+1 that maximizes the sum of instantenous utility and
    
74.     %   discounted continuation utility
    
75.     for h = 1 : gk
    
76.         for j = 1 : gz
    
77.             Tv(h,j) = max( u(h,:,j) + beta*(prh*v(:,1)' + prl*v(:,2)') );
    
78.         end
    
79.     end
    
80.     iter = iter + 1;
    
81.     diff = norm(Tv - v);
    
82.     v = Tv;
    
83. end
    
84. 
    
85. % Find the implicit decision rule for k_t+1 as a function of the state
    
86. %   variables k_t and z_t
    
87. for h = 1 : gk
    
88.     for j = 1 : gz
    
89.         % Using the [ ] syntax for max does not only give the value, but
    
90.         %   also the element chosen.
    
91.         [Tv,gridpoint] = max(u(h,:,j) + beta*(prh*v(:,1)' + prl*v(:,2)'));
    
92.         % Find what grid point of the k vector which is the optimal decision
    
93.         kgridrule(h,j) = gridpoint;
    
94.         % Find what value for k_t+1 which is the optimal decision
    
95.         kdecrule(h,j) = k(gridpoint);
    
96.     end
    
97. end
    
98. 
    
99. % Plot it and save it as a jpg-file
    
100. figure
     
101. plot(k,kdecrule,k,k);
     
102. xlabel('k_t')
     
103. ylabel('k_{t+1}')
     
104. print -djpeg kdecrule.jpg
     
105. 
     
106. % Compute the optimal decision rule for c as a function of the state
     
107. %    variables
     
108. for h = 1 : gk       % counter for k_t
     
109.     for j = 1 : gz   % counter for z_t
     
110.         cdecrule(h,j) = z(j)*k(h)^alpha + (1-delta)*k(h) - kdecrule(h,j);
     
111.     end
     
112. end
     
113. figure
     
114. plot(k,cdecrule)
     
115. xlabel('k_t')
     
116. ylabel('c_t')
     
117. print -djpeg decrulec.jpg
     
118. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
119. 
     
120. 
     
121. % Simulation
     
122. % Aribrary starting point
     
123. enostate = round(gk/2);
     
124. ksim(1) = k(enostate);
     
125. % Draw a sequence of 100 random variables and use the optimal decision rule
     
126. for i = 1 : 100
     
127.     draw = rand;
     
128.     if draw < prh
     
129.         exostate = 1;
     
130.     else
     
131.         exostate = 2;
     
132.     end
     
133.    
     
134.     kprimegrid = kgridrule(enostate,exostate);
     
135.     kprime = k(kprimegrid);
     
136.    
     
137.     zsim(i) = z(exostate);
     
138.     ksim(i+1) = kprime;
     
139.     ysim(i) = z(exostate)*ksim(i)^alpha;
     
140.     isim(i) = ksim(i+1) - (1-delta)*ksim(i);
     
141.     csim(i) = ysim(i) - isim(i);
     
142. end
     
143. figure
     
144. plot(1:100,ysim/mean(ysim),1:100,csim/mean(csim),1:100,isim/mean(isim),1:100,ksim(1:100)/mean(ksim(1:100)))
     
145. legend('y','c','i','k')
     
146. title('Simulation 100 draws, devations from mean')
     
147. print -djpeg simulation.jpg

复制代码

2# qibbxxt 你能提供程序，真是太感谢了。不知道你能不能提供这个程序来源的文章。

怎么给币法？我可以告诉你

我也很想知道这个程序的内容，能否告知？

这应该是你对所使用的工具箱中某工具的功能了解不全所致

请检查你的toolbox 确定你所使用的工具有你所使用的功能