如何在 MATLAB 中使用并行处理优化蒙特卡洛模拟？

问题描述 投票：0回答：1

我正在尝试优化经济模型的蒙特卡罗模拟，涉及三个 MATLAB 文件：StartMonCarFinal.m、SDsysFinal.m 和 SDMasterFinal.m。以下是每个文件功能的简要概述：

  • StartMonCarFinal.m：这是运行多个蒙特卡罗模拟的主脚本。它通过循环重复调用 SDMasterFinal.m，每次迭代代表一次模拟。
  • SDsysFinal.m：SDMasterFinal.m 调用的函数，处理一些内部系统变量的更新。
  • SDMasterFinal.m：执行模拟的核心脚本，调用 SDsysFinal.m 并处理大量数据和计算。

为了提高仿真的运行速度，我尝试将StartMonCarFinal.m中的常规for循环替换为parfor来实现并行处理。然而，在这次尝试中我遇到了几个问题：

  1. MATLAB 指示 save 命令不能在 parfor 循环内使用。我不知道如何克服这个限制。
  2. 我还担心SDMasterFinal.m中大量使用全局变量，这可能会与parfor的并行性发生冲突。

如果有任何关于如何有效并行处理此类蒙特卡罗模拟的建议，尤其是如何处理 parfor 循环中的保存问题和全局变量问题，我将不胜感激。任何有关如何提高此类脚本的并行效率的建议都会非常有帮助。

感谢您的时间和帮助！

  • StartMonCarFinal.m
niter = 100; %number of random realizations for montecarlo

for mcar = 1:niter
    
    save mcar mcar niter
    mcar
    SDMasterFinal
    
    load mcar
    
    if mcar > 1
        load OutputFinal
    end
    
    
    wMC(mcar,:,:) = w;
    Z_MMC(mcar,:,:) = Z_M;
    Z_SMC(mcar,:,:) = Z_S;
    Z_M_MorMC(mcar,:,:) = Z_M_Mor;
    Z_S_MorMC(mcar,:,:) = Z_S_Mor;
    c_MMC(mcar,:,:) = c_M;
    c_SMC(mcar,:,:) = c_S;
    p_MMC(mcar,:,:) = p_M;
    p_SMC(mcar,:,:) = p_S;
    MMC(mcar,:,:) = M;
    SMC(mcar,:,:) = S;
    L_MMC(mcar,:,:) = L_M;
    L_SMC(mcar,:,:) = L_S;
    L_MNANMC(mcar,:,:) = L_MNAN;
    L_SNANMC(mcar,:,:) = L_SNAN;
    LMC(mcar,:,:) = L;
    R_MMC(mcar,:,:) = R_M;
    R_SMC(mcar,:,:) = R_S;
    phi_MMC(mcar,:,:) = phi_M;
    phi_SMC(mcar,:,:) = phi_S;
    thetaMC(mcar,:,:) = theta;
    Po_MMC(mcar,:,:) = Po_M;
    Po_SMC(mcar,:,:) = Po_S;
    taoMC(mcar,:,:) = tao;
    tao_SMC(mcar,:,:) = tao_S;
    tao_MMC(mcar,:,:) = tao_M;
    ES_MMC(mcar,:,:) = ES_M;
    UbarMC(mcar,:,:) = Ubar;
    RIMC(mcar,:,:) = RI;
    Z_MvecMC(mcar,:,:) = Z_Mvec;
    Z_SvecMC(mcar,:,:) = Z_Svec;
    
    save OutputFinal *MC niter
    
end

save OutputFinal
  • SDMasterFinal.m

clear all 

global t h_M h_S alfa gamma mu sigma l a ksi1 ksi2 Lbar kappa
global w Z_M Z_S Z_M_Mor Z_S_Mor c_M c_S p_M p_S M S L_M L_S L R_M R_S phi_M phi_S theta Po_M Po_S tao ES_M Ubar RI

%Parameter Values
beta = .9747;
h_M = 0.9;
h_S = 1.4;
alfa = -1.5;
mu = .6;
gamma = 1-mu; %(needs to be 1-mu);
sigma =.6;
l = 1:500;
a = 45; %Pareto coefficient
ksi1 = 1.142793;
ksi2 = 0.0082433;
delta = 7.5;
Lbar = 100;
kappa = 0.08;
T = 240;

%Allocate variables
w = zeros(length(l),T);
Z_M = zeros(length(l),T);
Z_S = zeros(length(l),T);
Z_M_Mor = zeros(length(l),T);
Z_S_Mor = zeros(length(l),T);
c_M = zeros(length(l),T);
c_S = zeros(length(l),T);
p_M = zeros(length(l),T);
p_S = zeros(length(l),T);
M = zeros(length(l),T);
S = zeros(length(l),T);
L_M = zeros(length(l),T);
L_S = zeros(length(l),T);
L_MNAN = zeros(length(l),T);
L_SNAN = zeros(length(l),T);
L = zeros(length(l),T);
R_M = zeros(length(l),T);
R_S = zeros(length(l),T);
phi_M = zeros(length(l),T);
phi_S = zeros(length(l),T);
theta = zeros(length(l),T);
Po_M = zeros(length(l),T);
Po_S = zeros(length(l),T);
tao = zeros(length(l),T);
tao_S = zeros(length(l),T);
tao_M = zeros(length(l),T);
ES_M = zeros(length(l),T);
Ubar = zeros(1,T);
RI = zeros(1,T);
Z_Mvec = zeros(length(l),length(l));
Z_Svec = zeros(length(l),length(l));

 


X0(1) = 0.0492; %Ubar(1) initial value for solver
X0(2) = 1.035; %p_M(1,1) initial value for solver
X0(3) = 0.054; %RI(1) initial value for solver

for t = 1:T
    
    %t
    
    %Initial Productivity
    if t == 1
        Z_M_Mor(:,1) = .8 + .4*l/length(l);
        Z_S_Mor(:,1) = 1*ones(length(l),1);
    end
        
    F=@(R)SDsysFinal(R);
    options=optimset('MaxFunEvals',10000000000,'TolFun',0.0001,'MaxIter',1000000,'Display','off');
    [X,Eval] = fsolve(F,X0,options);
    
    X0(1) = Ubar(t);
    X0(2) = p_M(1,t);
    X0(3) = RI(t);
    
    %Eval

    %Diffusion 
    for i = 1:length(l)
        for j = 1:length(l)
            Z_Mvec(i,j) = exp(-delta*abs(i-j)/length(l))*Z_M(j,t);
            Z_Svec(i,j) = exp(-delta*abs(i-j)/length(l))*Z_S(j,t);
        end
        Z_M(i,t+1) = max(Z_Mvec(i,:));
        Z_S(i,t+1) = max(Z_Svec(i,:));
        
        %Save Morning Productivity
        Z_M_Mor(:,t+1) = Z_M(:,t+1);
        Z_S_Mor(:,t+1) = Z_S(:,t+1);
    end
    
end
  • SDsysFinal.m
function D = SDsysFinal(X)

global t h_M h_S alfa gamma mu sigma l a ksi1 ksi2 Lbar kappa 
global w Z_M Z_S Z_M_Mor Z_S_Mor c_M c_S p_M p_S M S L_M L_S L R_M R_S phi_M phi_S theta Po_M Po_S ES_M Ubar RI

p_S(:,t) = ones(length(l),1);

Ubar(t) = X(1);
p_M(1,t) = X(2);
RI(t) = X(3);

Z_M(:,t) = Z_M_Mor(:,t);
Z_S(:,t) = Z_S_Mor(:,t);

for i = 1:length(l)

    %Wages and consumption
    
    w(i,t) = Ubar(t)*((((p_M(i,t)./h_M)/(p_S(i,t)./h_S))^(1/(alfa-1)))*p_M(i,t) + p_S(i,t))...
        ./(h_M*(((p_M(i,t)./h_M)/(p_S(i,t)./h_S))^(alfa/(alfa-1))) + h_S).^(1/alfa) - RI(t);
        
    c_S(i,t) = (w(i,t) + RI(t))./(p_M(i,t).*(((p_M(i,t)./h_M)/(p_S(i,t)./h_S))^(1/(alfa-1)))+p_S(i,t));
    c_M(i,t) = c_S(i,t).*(((p_M(i,t)./h_M)/(p_S(i,t)./h_S))^(1/(alfa-1)));
    
    %Choose phi
    phi_M(i,t) = 1-((Z_M(i,t)*p_M(i,t)*gamma)/(ksi2*(a-1)*((w(i,t)/(mu*p_M(i,t)))^((1-gamma)/gamma))))^(-1/2);
    phi_S(i,t) = 1-((Z_S(i,t)*p_S(i,t)*gamma)/(ksi2*(a-1)*((w(i,t)/(sigma*p_S(i,t)))^((1-gamma)/gamma))))^(-1/2);
    
    if phi_M(i,t) < 0
        phi_M(i,t) = 0;
    end
    if phi_M(i,t) > 1;
        phi_M(i,t) = 1;
    end
    if phi_S(i,t) < 0
        phi_S(i,t) = 0;
    end
    if phi_S(i,t) > 1;
        phi_S(i,t) = 1;
    end

    %Production and labor choices
    L_M(i,t) = (w(i,t)./(mu*p_M(i,t).*((Z_M(i,t)*(phi_M(i,t)/(a-1)+1)).^gamma))).^(1/(mu-1));
    L_S(i,t) = (w(i,t)./(sigma*p_S(i,t).*((Z_S(i,t)*(phi_S(i,t)/(a-1)+1)).^gamma))).^(1/(sigma-1));
    
    if (ksi1 + ksi2*(1/(1-phi_M(i,t))))*w(i,t) > (((Z_M(i,t)).^gamma)*p_M(i,t).*(L_M(i,t).^mu)/((a-1)^gamma))...
            *(((phi_M(i,t)+a-1)^gamma - (a-1)^gamma));
        phi_M(i,t) = 0;
        L_M(i,t) = (w(i,t)./(mu*p_M(i,t).*(Z_M(i,t).^gamma))).^(1/(mu-1));
    end
    if (ksi1 + ksi2*(1/(1-phi_S(i,t))))*w(i,t) > (((Z_S(i,t)).^gamma)*p_S(i,t).*(L_S(i,t).^sigma)/((a-1)^gamma))...
            *(((phi_S(i,t)+a-1)^gamma - (a-1)^gamma));
        phi_S(i,t) = 0;
        L_S(i,t) = (w(i,t)./(sigma*p_S(i,t).*(Z_S(i,t).^gamma))).^(1/(sigma-1));
    end    
    
    %Land Use and Land Rents
    R_M(i,t) = p_M(i,t)*((Z_M(i,t)*(phi_M(i,t)/(a-1)+1)).^gamma).*(L_M(i,t).^mu)-w(i,t).*L_M(i,t)- w(i,t)*(ksi1+ksi2*(1/(1-phi_M(i,t))));
    R_S(i,t) = p_S(i,t)*((Z_S(i,t)*(phi_S(i,t)/(a-1)+1)).^gamma).*(L_S(i,t).^sigma)-w(i,t).*L_S(i,t)- w(i,t)*(ksi1+ksi2*(1/(1-phi_S(i,t))));
    
    
    if R_M(i,t) > R_S(i,t)
        theta(i,t) = 1; 
        phi_S(i,t) = 0;
    else
        theta(i,t) = 0;
        phi_M(i,t) = 0;
    end
    
    %Expected Outputs
    M(i,t) = ((Z_M(i,t)*(phi_M(i,t)/(a-1)+1)).^gamma).*L_M(i,t).^mu;
    S(i,t) = ((Z_S(i,t)*(phi_S(i,t)/(a-1)+1)).^gamma).*L_S(i,t).^sigma;
    
    %Aggregate Labor
    L(i,t) = theta(i,t).*L_M(i,t)+(1-theta(i,t)).*L_S(i,t);
    
    %Construct Excess Supply
    if i > 1
        if phi_M(i,t) > 0
            ES_M(i,t) =  ES_M(i-1,t) + (theta(i,t)*(M(i,t)- (w(i,t)/p_M(i,t))*(ksi1+ksi2*(1/(1-phi_M(i,t))))) - c_M(i,t).*L(i,t) - kappa*(abs(ES_M(i-1,t))))/length(l);
        else
            ES_M(i,t) =  ES_M(i-1,t) + (theta(i,t)*M(i,t) - c_M(i,t).*L(i,t) - kappa*(abs(ES_M(i-1,t))))/length(l);
        end    
    else
        if phi_M(i,t) > 0
            ES_M(i,t) =  (theta(i,t)*(M(i,t)- (w(i,t)/p_M(i,t))*(ksi1+ksi2*(1/(1-phi_M(i,t))))) - c_M(i,t).*L(i,t))/length(l);
        else
            ES_M(i,t) =  (theta(i,t)*M(i,t) - c_M(i,t).*L(i,t))/length(l);
        end
    end
    
    %Construct prices for i+1
    if ES_M(i,t) > 0
        if i < length(l)
            p_M(i+1,t) = p_M(i,t)*exp(kappa/length(l));
        end
    else
        if i < length(l)
            p_M(i+1,t) = p_M(i,t)*exp(-kappa/length(l));
        end    
    end
    
    %Generate Poisson draws
    Po_M(i,t) = rand;
    Po_S(i,t) = rand;
    
    %Generate Pareto and update productivity after innovation
    if Po_M(i,t) > phi_M(i,t);
        Z_M(i,t) = Z_M(i,t);
    else
        z = gprnd(1/a,1/a,1);
        Z_M(i,t) = z*Z_M(i,t); 
    end
    if Po_S(i,t) > phi_S(i,t);
        Z_S(i,t) = Z_S(i,t);
    else
        z = gprnd(1/a,1/a,1);
        Z_S(i,t) = z*Z_S(i,t); 
    end
        
    %Actual Outputs
    M(i,t) = theta(i,t).*(Z_M(i,t).^gamma).*L_M(i,t).^mu;
    S(i,t) = (1-theta(i,t)).*(Z_S(i,t).^gamma).*L_S(i,t).^sigma;
        
end

%Equilibrium in the goods market 
if (Ubar(t) > 0) %&& sum((1-theta(:,t)).*L_S(:,t))/length(l) < 100);
    d1 = ES_M(length(l),t);
else
    d1 = 1000000;
end    

%Equilibrium condition in labor market

if (p_M(1,t)> 0)% && p_M(1,t) < pp)
    d2 = sum(L(:,t))/length(l) - Lbar; 
else
    d2 = 1000000;
end 

%Calculate Rental Income (RI)

if (RI(t) > 0 && min(w(:,t)) > 0)
    d3 = RI(t) - (sum(theta(:,t).*R_M(:,t)) + sum((1-theta(:,t)).*R_S(:,t)))/(Lbar*length(l));  
else
    d3 = 1000000;
end

D = [d1,d2,d3];
matlab montecarlo parfor
1个回答
0
投票

MATLAB 指示 save 命令不能在 parfor 循环内使用。我不知道如何克服这个限制。

您可以创建一个函数来在 forloop 内进行保存，如下所示：

function parsave(fname, x,y)
  save(fname, 'x', 'y')
end

我还担心SDMasterFinal.m中大量使用全局变量，这可能会与parfor的并行性发生冲突。

我建议根本不要使用全局，因为它很难调试，我建议尽可能编辑代码以使用函数，以便更容易实现 MC。

您可以创建一个名为 

states
的结构，其中包含所有所需的变量并将其包含在函数中。我建议使用较小的结构，以便您确切地知道函数内部和外部的内容

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.