程序是由卡内基梅隆大学的Dan Li 编写的，邮箱是danli@andrew.cmu.edu，个人网页：www.danli.org
此贴仅供学习参考。
假设没有人口增长和技术进步，劳动力的供给是没有弹性的。模型是Deterministic neoclassical growth model。
所用的生产函数是科布道格拉斯函数，效用函数是幂指数效用函数。
先定义生产函数：
% Cobb-Douglas production function
% Y=F(K,N)=K^{alpha}*N^{1-alpha}
% In per capita form y=k^{alpha}

function y=cobb(k)
global alpha
if alpha>0 & alpha<1
y=k.^alpha;
else
y=-inf*abs(k);
disp('invalid capital share')
end

再定义效用函数：
% Power utility function u(c)=c^(1-gamma)/(1-gamma) if gamma ~= 1;
% u(c)=ln(c) if gamma=1.

% Utility is only defined if c>=0 for gamma<1 and c>0 if gamma>1

% The risk aversion coefficient is pass on to the function by
% define it as a global variable.


function u=util(c)

global gamma;
u=-inf*ones(size(c)); % initialize into -inf, change only those c positive
cpos=find(c>0); % find those positive c

if gamma==1
u(cpos)=log(c(cpos));
else
u(cpos)=c(cpos).^(1-gamma)./(1-gamma);
end

然后定义画图函数：
function vgc_plot(k,kstar,v,Polc,g, METHOD)
% plot the value function against k
% plot capital policy for next period against captial this period
% save variable
% input is the k vector, optimal-- kstar, value --v, capital policy - g
% METHOD=1 if deterministic noninterpolation
% Method=2----Deterministic, continuous state space model, With interpolation
gstar=ppval(spline(k,g),kstar);% g value at optimal kstar
vstar=ppval(spline(k,v),kstar);% v value at optimal kstar

T=char( ' No Interpolation, Deterministic',...
' With Interpolation, Deterministic');

figure(1);
plot(k,v); hold on;
plot(kstar,vstar,'ro');
text(kstar+.1,vstar-.1,['k*=' num2str(kstar)]);
title(['Value Function--' T(METHOD,:)]);
xlabel('capital k');ylabel('value');
axis square; hold off;
print('-dpsc',['ValFun' num2str(METHOD)]);

figure(2);
subplot(2,1,1); plot(k,Polc);
title('Policy Function-Consumption');
xlabel('Capital Kt'); ylabel('Consumption');

subplot(2,1,2); plot(k,g,k,k,'-'); hold on;
plot(kstar,gstar,'ro');
text(kstar+.1,gstar-.1,['g*=' num2str(gstar) '| k*=' num2str(kstar)]);
title('Policy Function-Capital');
xlabel('Capital Kt'); ylabel('Capital Next Period');
hold off;

print('-dpsc', ['PolFun' num2str(METHOD)]);
save(['growth' num2str(METHOD)]);

最后运行新古典增长模型：
clc;
clear all;

% ------Model Parameter specification---------
global gamma alpha
gamma=1; % risk aversion coefficient, log utility function
alpha=0.3; % capital share
beta=0.96; % discount factor
delta=.1; % depreciation rate
z=1; % initial state for deteriministic case

method=input(['Method=?,\n'...
'1 for Deterministic,No interpolation \n' ...
'2 for Deterministic With Interpolation,\n']);

%------ Parameter of computation--------------------------
nk=input(['Number of grid?(20-200 for Method 1)\n'...
'(<=30 for Method 2, <=200 for Method 3)>>']); % number of grid
maxloop=40; % number of maximum iteration
tol=10^(-3); % tolerence level for value function iteration

tic; % start the clock

% ******Low and upper bound of state variable k ******
% Steady state value kstar, f'(kstar)=1/beta-1+delta

kstar=(((1/(alpha*beta))-((1-delta)/alpha)))^(1/(alpha-1));
lowk=0.8*kstar;
upk=1.2*kstar;

k=linspace(lowk,upk,nk)';% build up grid
[K,Kprime]=meshgrid(k);

v=zeros(size(k))*ones(size(z)); vn=v;
g=v; gn=g;
loop=0; dv=1; dk=1;

if method==1
%--------------------------------------
%No Interpolation Method, Method=1
%---------------------------------------
C=cobb(K)+(1-delta)*K-Kprime;
U=util(C);

kpos=100*ones(size(k));% initial value of loop and distance of value

while dv>tol & loop<=maxloop & dk
loop=loop+1;
[vn,kposn]=max(U+beta*v*ones(1,nk));%new value vector
vn=vn'; kposn=kposn';
dv=max(abs((vn-v)./(v+eps)));
dk=~all(kpos==kposn); % whether or not position exactly the same
v=vn; kpos=kposn;
end
g=k(kpos);%policyfunction of next period capital

elseif method==2;
%--------------------------------------
%Interpolation Method, Method=2
%---------------------------------------
while dv>tol & dk>tol & loop<=maxloop
loop=loop+1;
pp = spline(k,v);
vf=inline('-(util(-y+cobb(k)+(1-delta)*k)+beta*ppval(pp,y))','y','k','delta','beta','pp');
for ik=1:nk
[gn(ik),vn(ik)]=fminbnd(vf,lowk,upk,[],k(ik),delta,beta,pp) ;
vn(ik)=-vn(ik);
end;
dk=max(abs((gn-g)./(g+eps)));
dv=max(abs((vn-v)./(v+eps)));
v=vn; g=gn;
end;
end
%------------------------------------------------------------
% use calculated k g to get consumption policy and plot , save
%--------------------------------------------------------------
%Polc=cobb(k)*z+(1-delta)*k*ones(size(z))-g; %policy function of consumption
Polc=cobb(k)+(1-delta)*k-g; %policy function of consumption

vgc_plot(k,kstar,v,Polc,g, method);

%--------------so other output for observation of convergence--------------
t=toc% read clock