%------------ THE STOCHASTIC MODEL: DUOPOLY and SOCIAL PLANNER----------------%
% Date: October. 10, 2010 %
 
% ATTENTION: 
% 1. with "con" we indicate the PUBLIC SUBSIDY.
% 2. with "theta" we indicate A in the text.
 
clear;
clc;
 
M = 10 ;     % Number of simulations
N=50;        % Number of simulation periods (time)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. INITIALIZATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.1 Initialize of the simulations (Matrices for the variables)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% Subsidy
con_sim_MC = zeros(N-1,M);  % Subsidy initialization for the matrix containing the M simulations
con_sim_sp_MC = zeros(N-1,M);
 
% Investment in R&D
inv2_sim_MC = zeros(N-1,M);
inv2_sim_sp_MC = zeros(N-1,M);
 
% R&D Capital
k_sim_MC = zeros(N,M);
k_sim_sp_MC = zeros(N,M);
 
% Ratio S/R
rapp_MC = zeros(N-1,M);
rapp_sp_MC = zeros(N-1,M);
 
% Profit Agency
prof_ag_MC = zeros(N-1,M);
prof_ag_sp_MC = zeros(N-1,M);
 
% Profit Firm
prof_imp_MC = zeros(N-1,M);
prof_imp_sp_MC = zeros(N-1,M);
 
 
for index=1:M;                                      
 
%------------------------------------------------------------
% 1.2 Set parameters for the Bellman equation (see the paper)
%------------------------------------------------------------
 
delta=0.15; % Depreciation rate
beta0=0.96; % Discount factor
beta = 8; 
gamma= 1;
mu=1.5
fi0=5;
fi1=25;
xi=(1-delta)/(2*fi1*(1-delta)+2*gamma);
fi=(1-delta-xi*gamma)/(1-delta);
a= (fi0-mu)/beta; 
b= (1-delta)/(beta*xi); % Expressed in absolute value
c=fi/(beta*xi); % Expressed in absolute value
psi =20; % Agency parameter
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.3 Set parameters for the first-order Markov process of A=theta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
rho=0.5; % Parameter in first-order Markov Process --> IT INDICATES STATE PERSISTENCY
epsilon = 1         % Stochastic shock
thetaA = + epsilon  % Positive schock
thetaB = - epsilon  % Negative schock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. VALUE FUNCTION ITERATION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%----------------------------------------
% 2.1 Set State-Space 
%----------------------------------------
 
n=100; % Number of grids
kmax=1; % Max k. 
% Mathematically it is defined as:
k=(kmax/n:kmax/n:kmax)' ; % State space: sequence of n possible states
i=ones(n,1);
k2=(i*k');
[I1,J1]=find(k2'-(1-delta)*k2<0);  % Constraint for investment in R&D >=0
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[I2a,J2a]=find(-a/thetaA*ones(n,n)-b/thetaA*k2+c/thetaA*k2'<0);   % Constraint for subsidy >=0
Ia=[I1;I2a];                                                      % Vertical concatenation
Ja=[J1;J2a];
[I2b,J2b]=find(-a/thetaB*ones(n,n)-b/thetaB*k2+c/thetaB*k2'<0);
Ib=[I1;I2b];                                                      % Vertical concatenation
Jb=[J1;J2b];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the "DUOPOLY (COMPETITIVE EQUILIBIUM):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.2 Calculate the “policy function” 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%-------------------------------------------------
% Value function Iteration of the Bellman Equation 
%-------------------------------------------------
 
j=0; % Index for the loop to keep track of number of iterations
h=1; % Index [parameter(arbitrary number used for convergence)]
vA=zeros(n,1);
vB=zeros(n,1);
while h > 0.005; % Criterion for convergence
    k3A=((-a/thetaA*ones(n,n)-b/thetaA*k2+c/thetaA*k2').*(ones(n,n)+k2'-(1-delta)*k2))-psi*((-a/thetaA*ones(n,n)-b/thetaA*k2+c/thetaA*k2').^(2))+beta0*(vA*i'*((1+rho)/2)+vB*i'*((1-rho)/2)); %k3A is a matrix
    k3B=((-a/thetaB*ones(n,n)-b/thetaB*k2+c/thetaB*k2').*(ones(n,n)+k2'-(1-delta)*k2))-psi*((-a/thetaB*ones(n,n)-b/thetaB*k2+c/thetaB*k2').^(2))+beta0*(vA*i'*((1+rho)/2)+vB*i'*((1-rho)/2)); %k3B is a matrix
 
    for l=1:length(Ja);
        k3A( Ia(l),Ja(l) )=-10000;  % Penalizing the grids where c & i are negative                              
    end;
    for l=1:length(Jb);
        k3B( Ib(l),Jb(l) )=-10000;  % Penalizing the grids where c & i are negative                              
    end;
        
    [maxvecA,locA]=max(k3A);   % Max value per column; maxvec=values; loc=location (row)
    [maxvecB,locB]=max(k3B);   % Max value per column; maxvec=values; loc=location (row)
 
    h=max(max(abs(maxvecA'-vA)), max(abs(maxvecB'-vB))) ; 
    vA=maxvecA';
    vB=maxvecB';
    j=j+1;
end;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.3 Generate the “policy function” and graph it
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
polA=zeros(n,1);        % Initializing the policy function 
for q=1:n;
    polA(q)=k(locA(q));  % Generating the policy function
end;
 
polB=zeros(n,1);        % Initializing the policy function 
for q=1:n;
    polB(q)=k(locB(q));  % Generating the policy function
end;
 
% figure (1);
% %% VALUE FUNCTION
% figure(1);
% plot (maxvecA, 'DisplayName', 'maxvecA', 'YDataSource', 'maxvecA'); hold all; plot (maxvecB, 'DisplayName', 'maxvecB', 'YDataSource', 'maxvecB'); hold off;
% legend('V(thetaA)', 'V(thetaA)');
% title('VALUE FUNCTION');
% xlabel('CURRENT CAPITAL');
% ylabel('VALUE FUNCTION');
 
%% POLICY FUNCTION
% figure(2);
% plot (polA, 'DisplayName', 'polA', 'YDataSource', 'polA'); hold all; plot (polB, 'DisplayName', 'polB', 'YDataSource', 'polB'); hold off; 
% legend('Policy(thetaA)', 'Policy(thetaB)');
% title('POLICY FUNCTION');
% xlabel('CAPITAL');
% ylabel('POLICY FUNCTION');
 
%% INVESTMENT (R&D)
% figure(4);
% invA=polA-(1-delta)*k;
% invB=polB-(1-delta)*k;
% plot (invA, 'DisplayName', 'invA', 'YDataSource', 'invA'); hold all; plot (invB, 'DisplayName', 'invB', 'YDataSource', 'invB'); hold off; 
% legend('Investment(thetaA)', 'Investment(thetaB)');
% ylabel('Investment');
% xlabel('Capital');
% title('INVESTMENT');
 
%% SUBSIDY 
% figure(5);
% conA= -a/thetaA-b/thetaA*k+c/thetaA*polA; 
% conB= -a/thetaB-b/thetaB*k+c/thetaB*polB; 
% plot (conA, 'DisplayName', 'conA', 'YDataSource', 'conA'); hold all; plot (conB, 'DisplayName', 'conB', 'YDataSource', 'conB'); hold off;
% legend('Subsidy(thetaA)', 'Subsidy(thetaB)');
% ylabel('Subsidy');
% xlabel('Capital');
% title('SUBSIDY');
 
 



%---------------------
% 3. SIMULATION
%---------------------
 
%---------------------------------------------
% 3.1 Simulate the Markov Chain of A=theta
%---------------------------------------------
% epsilon=1;   % Parameter defining the value of theta
% thetaA= - epsilon;
% thetaB= + epsilon;
% rho=0.5
 
P = [(1+rho)/2,(1-rho)/2;
     (1-rho)/2,(1+rho)/2]; %% n-by-n "transition matrix"
 
pi0 = [1,0];     %% Initial probability distribution: we start with "tethaA=positive-shock"
 
x   = [thetaA;thetaB];     %% n-by-1 state
 
%%%%% Markov-chain Simulation:
 
%[chain,states] = simulate_markov(x,P,pi0,T);
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.2 Set the vector for the simulated time-series of the variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%N=50;  
k_sim=zeros(N,1);      % Vector of simulated time series of capital
inv_sim=zeros(N-1,1);
con_sim=zeros(N-1,1);
inv2_sim=zeros(N-1,1);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.3 Choose initial values for the R&D capital and then simulate 
%     the temporal movement of all the variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
L=10;            % Location of initial level of current capital where simulation starts
k_sim(1)=k(L);   % Note: subsidy and investment in R&D cannot be computed yet
 
[chain,states] = simulate_markov2(x,P,pi0,N);
%figure(1)
%plot(chain,'bo-')
%title('Markov chain')
 
for i=2:N;
    if chain(i)==thetaA;
    k_sim(i) = polA(L); 
    L= find(k==polA(L)); 
    inv_sim(i-1)=k_sim(i)-(1-delta)*k_sim(i-1);
    con_sim(i-1)=-a/thetaA-b/thetaA*k_sim(i-1)+c/thetaA*k_sim(i);  % subsidy
    inv2_sim(i-1)=xi*(fi0-mu)+beta*xi*thetaA*con_sim(i-1)+xi*gamma/(1-delta)*k_sim(i)
    elseif chain(i)==thetaB;
    k_sim(i) = polB(L); 
    L= find(k==polB(L));
    inv_sim(i-1)=k_sim(i)-(1-delta)*k_sim(i-1);
    con_sim(i-1)=-a/thetaB-b/thetaB*k_sim(i-1)+c/thetaB*k_sim(i);  % subsidy
    inv2_sim(i-1)=xi*(fi0-mu)+beta*xi*thetaB*con_sim(i-1)+xi*gamma/(1-delta)*k_sim(i)
    end
end
 
% Ratio between subsidy and R&D investment
rapp = con_sim./inv2_sim 
 
% Profit Firm
prof_imp = inv2_sim.*(fi0*ones(N-1,1)-fi1*inv2_sim)-(mu*ones(N-1,1)-beta*con_sim-gamma*k_sim(2:N)).*inv2_sim 
 
% Profit Agency
prof_ag = con_sim.*(ones(N-1,1)+inv2_sim)-psi*(con_sim).^(2)
 
% figure(5);
% plot (k_sim);
% Ylabel('capital');
% xlabel('simulation');
% title('Simulated R&D CAPITAL STOCK function');
% 
% figure(6);
% plot(con_sim);
% ylabel('Subsidy');
% title('Simulated SUBSIDY function');
% 
% figure(7);
% plot(prof_imp);
% ylabel('Firm profit');
% title('Simulated FIRM PROFIT function');
% 
% figure(8);
% plot(prof_ag);
% ylabel('Agency profit');
% title('Simulated AGENCY PROFIT function');
% 
% figure(9);
% plot(rapp);
% ylabel('Ratio S/R');
% title('Simulated RATIO of S to R function');
% 
% %figure(10);
% %plot(inv_sim);
% %ylabel('investment');
% %title('Simulated investment function');
% 
% figure(11);
% plot(inv2_sim);
% ylabel('Firm R&D investment');
% title('Simulated FIRM R&D INVESTMENT function (2)');
% 
% figure(12)
% plot(chain,'bo-')
% title('Simulated Markov chain')
 
% End Duopoly

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the SOCIAL PLANNER (COOPERATIVE EQULIBRIUM):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. VALUE FUNCTION ITERATION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%----------------------------------------
% 2.1 Set State-Space 
%----------------------------------------
 
n=100;  % Number of grids
kmax=1; % Max k. 
% Mathematically it is defined as:
k=(kmax/n:kmax/n:kmax)' ; % State space: sequence of n possible states
i=ones(n,1);
k2=(i*k');
[I1,J1]=find(k2'-(1-delta)*k2<0);  % Constraint for investment in R&D >=0
 
 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.2 Calculate the “policy function” 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%-------------------------------------------------
% Value function Iteration of the Bellman Equation 
%-------------------------------------------------
 
j=0; % Index for the loop to keep track of number of iterations
h=1; % Index [parameter (arbitrary number used for convergence)]
vA=zeros(n,1);
vB=zeros(n,1);
while h > 0.005; % Criterion for convergence
    
    k3A=((-a/thetaA*ones(n,n)-b/thetaA*k2+c/thetaA*k2').*(ones(n,n)+k2'-(1-delta)*k2))-psi*((-a/thetaA*ones(n,n)-b/thetaA*k2+c/thetaA*k2').^(2)) ...
    + (k2'-(1-delta)*k2).*(fi0*ones(n,n)-fi1*(k2'-(1-delta)*k2))-(mu*ones(n,n)-beta*thetaA*(-a/thetaA*ones(n,n)-b/thetaA*k2+c/thetaA*k2')-gamma*k2).*(k2'-(1-delta)*k2)+beta0*(vA*i'*((1+rho)/2)+vB*i'*((1-rho)/2));
 
    k3B=((-a/thetaB*ones(n,n)-b/thetaB*k2+c/thetaB*k2').*(ones(n,n)+k2'-(1-delta)*k2))-psi*((-a/thetaB*ones(n,n)-b/thetaB*k2+c/thetaB*k2').^(2)) ...
    + (k2'-(1-delta)*k2).*(fi0*ones(n,n)-fi1*(k2'-(1-delta)*k2))-(mu*ones(n,n)-beta*thetaB*(-a/thetaB*ones(n,n)-b/thetaB*k2+c/thetaB*k2')-gamma*k2).*(k2'-(1-delta)*k2)+beta0*(vA*i'*((1+rho)/2)+vB*i'*((1-rho)/2));
   
    for l=1:length(Ja);
        k3A( Ia(l),Ja(l) )=-10000;    % Penalizing the grids where c & i are negative                              
    end;
    for l=1:length(Jb);
        k3B( Ib(l),Jb(l) )=-10000;    % Penalizing the grids where c & i are negative                              
    end;
        
    [maxvecA,locA]=max(k3A);   % Max value per column; maxvec=values; loc=location (row)
    [maxvecB,locB]=max(k3B);   % Max value per column; maxvec=values; loc=location (row)
 
    h=max(max(abs(maxvecA'-vA)), max(abs(maxvecB'-vB))) ; 
    vA=maxvecA';
    vB=maxvecB';
    j=j+1;
end;
 
%--------------------------------------------------------
% 2.3 Generate the “policy function” and graph it
%--------------------------------------------------------
 
polA_sp=zeros(n,1);        % Initializing the policy function 
for q=1:n;
    polA_sp(q)=k(locA(q));  % Generating the policy function
end;
 
polB_sp=zeros(n,1);        % Initializing the policy function 
for q=1:n;
    polB_sp(q)=k(locB(q));  % Generating the policy function
end;
 
% figure (1);
% %% VALUE FUNCTION
% figure(1);
% plot (maxvecA, 'DisplayName', 'maxvecA', 'YDataSource', 'maxvecA'); hold all; plot (maxvecB, 'DisplayName', 'maxvecB', 'YDataSource', 'maxvecB'); hold off;
% legend('V(thetaA)', 'V(thetaA)');
% title('VALUE FUNCTION');
% xlabel('CURRENT CAPITAL');
% ylabel('VALUE FUNCTION');
 
% %% POLICY FUNCTION
% figure(2);
% plot (polA_sp, 'DisplayName', 'polA', 'YDataSource', 'polA'); hold all; plot (polB_sp, 'DisplayName', 'polB', 'YDataSource', 'polB'); hold off; 
% legend('Policy(thetaA)', 'Policy(thetaB)');
% title('POLICY FUNCTION');
% xlabel('CAPITAL');
% ylabel('POLICY FUNCTION');
 
%% INVESTMENT (R&D)
% figure(4);
% invA_sp=polA_sp-(1-delta)*k;
% invB_sp=polB_sp-(1-delta)*k;
% plot (invA_sp, 'DisplayName', 'invA', 'YDataSource', 'invA'); hold all; plot (invB_sp, 'DisplayName', 'invB', 'YDataSource', 'invB'); hold off; 
% legend('Investment(thetaA)', 'Investment(thetaB)');
% ylabel('Investment');
% xlabel('Capital');
% title('INVESTMENT');
 
% %% SUBSIDY 
% figure(5);
% conA_sp= -a/thetaA-b/thetaA*k+c/thetaA*polA_sp; 
% conB_sp= -a/thetaB-b/thetaB*k+c/thetaB*polB_sp; 
% plot (conA_sp, 'DisplayName', 'conA', 'YDataSource', 'conA'); hold all; plot (conB_sp, 'DisplayName', 'conB', 'YDataSource', 'conB'); hold off;
% legend('Subsidy(thetaA)', 'Subsidy(thetaB)');
% ylabel('Subsidy');
% xlabel('Capital');
% title('SUBSIDY');
 
%---------------------
% 3. SIMULATION
%---------------------
 
%---------------------------------------------
% 3.1 Simulate the Markov Chain of A=theta
%---------------------------------------------
 
%epsilon=0.1;       % Parameter defining the value of theta
%thetaA=1-epsilon;
%thetaB=1+epsilon;
%rho=0.9
 
%P = [(1+rho)/2,(1-rho)/2;
%     (1-rho)/2,(1+rho)/2]; %% n-by-n "transition matrix"
 
%pi0 = [1,0];     %% Initial probability distribution
 
%x   = [thetaA;thetaB];     %% n-by-1 state
 
%%%%% Markov-chain Simulation:
 
%[chain,states] = simulate_markov(x,P,pi0,T);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.2 Set the vector for the simulated time-series of the variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% N=50;
k_sim_sp=zeros(N,1);      % Vector of simulated time series of capital
inv_sim_sp=zeros(N-1,1);
con_sim_sp=zeros(N-1,1);
inv2_sim_sp=zeros(N-1,1);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.3 Choose initial values for the R&D capital and then simulate 
%     the temporal movement of all the variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
L=10;               % Location of initial level of current capital where simulation starts
k_sim_sp(1)=k(L);   % Note: subsidy and investment cannot be computed yet
 
%[chain,states] = simulate_markov2(x,P,pi0,N);
%figure(1)
%plot(chain,'bo-')
%title('Markov chain')
 
for i=2:N;
    if chain(i)==thetaA;
    k_sim_sp(i) = polA_sp(L); 
    L= find(k==polA_sp(L)); 
    inv_sim_sp(i-1)=k_sim_sp(i)-(1-delta)*k_sim_sp(i-1);
    con_sim_sp(i-1)=-a/thetaA-b/thetaA*k_sim_sp(i-1)+c/thetaA*k_sim_sp(i);  % subsidy
    inv2_sim_sp(i-1)=xi*(fi0-mu)+beta*xi*thetaA*con_sim_sp(i-1)+xi*gamma/(1-delta)*k_sim_sp(i)
    elseif chain(i)==thetaB;
    k_sim_sp(i) = polB_sp(L); 
    L= find(k==polB_sp(L));
    inv_sim_sp(i-1)=k_sim_sp(i)-(1-delta)*k_sim_sp(i-1);
    con_sim_sp(i-1)=-a/thetaB-b/thetaB*k_sim_sp(i-1)+c/thetaB*k_sim_sp(i);  % subsidy
    inv2_sim_sp(i-1)=xi*(fi0-mu)+beta*xi*thetaB*con_sim_sp(i-1)+xi*gamma/(1-delta)*k_sim_sp(i)
    end
end
 
% Ratio between subsidy and R&D investment
rapp_sp = con_sim_sp./inv2_sim_sp 
 
% Profit Firm
prof_imp_sp = inv2_sim_sp.*(fi0*ones(N-1,1)-fi1*inv2_sim_sp)-(mu*ones(N-1,1)-beta*con_sim_sp-gamma*k_sim_sp(2:N)).*inv2_sim_sp 
 
% Profitto Agency
prof_ag_sp = con_sim_sp.*(ones(N-1,1)+inv2_sim_sp)-psi*(con_sim_sp).^(2)
 
% figure(5);
% plot (k_sim_sp);
% Ylabel('capital');
% xlabel('simulation');
% title('Simulated R&D CAPITAL STOCK function SOCIAL PLANNER');
% 
% figure(6);
% plot(con_sim_sp);
% ylabel('Subsidy');
% title('Simulated SUBSIDY function SOCIAL PLANNER');
% 
% figure(7);
% plot(prof_imp_sp);
% ylabel('Firm profit');
% title('Simulated FIRM PROFIT function SOCIAL PLANNER');
% 
% figure(8);
% plot(prof_ag_sp);
% ylabel('Agency profit');
% title('Simulated AGENCY PROFIT function SOCIAL PLANNER');
% 
% figure(9);
% plot(rapp_sp);
% ylabel('Ratio S/R');
% title('Simulated RATIO of S to R function SOCIAL PLANNER');
% 
% %figure(10);
% %plot(inv_sim_sp);
% %ylabel('investment');
% %title('Simulated investment function SOCIAL PLANNER');
% 
% figure(11);
% plot(inv2_sim_sp);
% ylabel('Firm R&D investment');
% title('Simulated FIRM R&D INVESTMENT function (2) SOCIAL PLANNER');
% 
% figure(12)
% plot(chain,'bo-')
% title('Simulated Markov chain')
% End
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Overlapped graphs for DUOPOLY and SOCIAL-PLANNER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% figure(13)
% plot (inv2_sim,':','LineWidth',1.5); hold all; plot (inv2_sim_sp, 'LineWidth',1.5); hold off; 
% ylabel('R&D investment');
% xlabel('Time')
% title('Simulated R&D INVESTMENT');
% legend('Duopoly' , 'Social-planner')
% 
% 
% figure(14)
% plot (con_sim,':','LineWidth',1.5); hold all; plot (con_sim_sp, 'LineWidth',1.5); hold off; 
% ylabel('Subsidy');
% xlabel('Time')
% title('Simulated SUBSIDY');
% legend('Duopoly' , 'Social-planner')
% 
% 
% figure(15)
% plot (rapp,':','LineWidth',1.5); hold all; plot (rapp_sp, 'LineWidth',1.5); hold off; 
% ylabel('S/R');
% xlabel('Time')
% title('Simulated RATIO of S to R');
% legend('Duopoly' , 'Social-planner')
% 
% 
% figure(16)
% plot (k_sim,':','LineWidth',1.5); hold all; plot (k_sim_sp, 'LineWidth',1.5); hold off; 
% ylabel('R&D capital');
% xlabel('Time')
% title('Simulated R&D CAPITAL');
% legend('Duopoly' , 'Social-planner')
% 
% 
% figure(17)
% plot (prof_imp,':','LineWidth',1.5); hold all; plot (prof_imp_sp, 'LineWidth',1.5); hold off; 
% ylabel('Firm profit');
% xlabel('Time')
% title('Simulated FIRM PROFIT');
% legend('Duopoly' , 'Social-planner')
% 
% 
% figure(18)
% plot (prof_ag,':','LineWidth',1.5); hold all; plot (prof_ag_sp, 'LineWidth',1.5); hold off; 
% ylabel('Agency profit');
% xlabel('Time')
% title('Simulated AGENCY PROFIT');
% legend('Duopoly' , 'Social-planner')
 
% End
 
%---------------------------------------
% CALCULATION (NxM MATRICES)
%---------------------------------------
 
% Subsidy
con_sim_MC(:,index) = con_sim;
con_sim_sp_MC(:,index) = con_sim_sp;
 
% R&D investment
inv2_sim_MC(:,index) = inv2_sim;
inv2_sim_sp_MC(:,index) = inv2_sim_sp; 
 
% R&D Capital
k_sim_MC(:,index) = k_sim;
k_sim_sp_MC(:,index) = k_sim_sp;
 
% Ratio S/R
rapp_MC(:,index) = rapp;
rapp_sp_MC(:,index) = rapp_sp;
 
% Profit Agency
prof_ag_MC(:,index) = prof_ag;
prof_ag_sp_MC(:,index) = prof_ag_sp;
 
% Profit Firm
prof_imp_MC(:,index) = prof_imp;
prof_imp_sp_MC(:,index) = prof_imp_sp;
 
end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.4 Calculate variables' average over the M Monte-Carlo simulations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% Subsidy
con_sim_MC_mean = mean(con_sim_MC,2);
con_sim_sp_MC_mean = mean(con_sim_sp,2);
 
% R&D investment
inv2_sim_MC_mean = mean(inv2_sim_MC,2);
inv2_sim_sp_MC_mean = mean(inv2_sim_sp_MC,2); 
 
% R&D Capital
k_sim_MC_mean = mean(k_sim_MC,2);
k_sim_sp_MC_mean = mean(k_sim_sp_MC,2);
 
% Ratio S/R
rapp_MC_mean = mean(rapp_MC,2);
rapp_sp_MC_mean = mean(rapp_sp_MC,2);
 
% Profit Agency
prof_ag_MC_mean = mean(prof_ag_MC,2);
prof_ag_sp_MC_mean = mean(prof_ag_sp_MC,2);
 
% Profit Firm
prof_imp_MC_mean = mean(prof_imp_MC,2);
prof_imp_sp_MC_mean = mean(prof_imp_sp_MC,2);
 
% Welfare 
welfare_MC_mean = prof_ag_MC_mean + prof_imp_MC_mean;
welfare_sp_MC_mean = prof_ag_sp_MC_mean + prof_imp_sp_MC_mean;
 
% Marginal rate of return (MRR)
mrr_MC_mean = fi0 - fi1*inv2_sim_MC_mean;
mrr_sp_MC_mean = fi0 - fi1*inv2_sim_sp_MC_mean;
 
% Marginal capital costs (MCC)
mcc_MC_mean = mrr_MC_mean - (prof_imp_MC_mean./inv2_sim_MC_mean);
mcc_sp_MC_mean = mrr_sp_MC_mean - (prof_imp_sp_MC_mean./inv2_sim_sp_MC_mean);
 
% END OF SIMULATION
 
% Percentage calculation of the difference between the "ratio S/R"
% for the Social-planner and the "ratio S/R" in the Duopoly. 
 
scarto_rapp =  ((rapp_sp_MC_mean - rapp_MC_mean)./rapp_MC_mean)*100
mean(scarto_rapp)
 
% END
