clear all
clc

%*This program simulates datasets following DGP1 in ACF (2015)*/

rng(1000);

nreps = 1000; %* Number of MC replications */

%* To create DGP1 */

siglnw = 0.1; %* Measures SD of firm specific wages - For the original DGP1, this should be set to 0.1. We also modify the DGP1 by increasing the standard deviation of the log wage to 0.5 */
timeb = 0.5; %/* Determines Point in time when labor is chosen*/
sigoptl = 0; %/* Optimization error in labor. For DGP1, this should be set to zero */

r = 1;
while r<=nreps;   %/* Loop over the MC replications */
r

%/* STEP 1:  GENERATE THE DATA */


n = 1000;   %/* Number of firms */

%/* We decided to take data from an approximate steady state of the dynamic model.  We start firms with (close to) zero level of capital, and run the model for 100 periods (which seemed to get us pretty close to the steady state).  The data used in estimation is the last 10 of these periods */

overallt = 100;  %/* Number of time periods in simulation in total (~100 seemed to get to approximate steady state) */
starttime =90;   %/* Time period at which observed data starts */
t = 10;          %/* Time periods in observed data */


%/* Matrices to store data */

lnkdata = zeros(n,overallt);   %/* ln(capital) */
lnldata = zeros(n,overallt);   %/* ln(labor)   */
lnmdata = zeros(n,overallt);   %/* ln(intermediate input) */
lnwdata = zeros(n,overallt);   %/* ln(wage) */
lnpdata = zeros(n,overallt);   %/* ln(output price) */
lnydata = zeros(n,overallt);   %/* ln(output) */
omgdata = zeros(n,overallt);  %/* omega(t) */
omgdataminusb = zeros(n,overallt); %/* omega(t-b) */
epsdata = zeros(n,overallt);  %/* epsilon */


%/* Production Function Parameters */

alpha0 = 0; %log(alpha0) = log(1) = 0
alphal = 0.6;
alphak = 0.4;

%/* Generate epsilons */

sigeps = 0.1;
epsdata = sigeps.*randn(n,overallt); %epsdata = sigeps*rndn(n,overallt); 

%/* Generate omega and ln(wage) processes */

sigomg = 0.3;  %/* SD of omega */

rho = 0.7; %/* AR(1) coefficient for omega */

rholnw = 0.3; %/* AR(1) coefficient for ln(wage) - note: SD defined at top of program */

%/* The next five lines define parameters necessary to subdivide the omega process into omega(t-1) -> omega(t-b) -> omega(t).  Note that this depends on the timeb (b) parameter, defined at top of program */

rhofirst = rho^(1-timeb);
rhosecond = rho^timeb;

sigxi = sqrt((1-rho^2)*sigomg^2);  %/* SD of innovation in omega - set such that variance of omega(t) is constant across time */
sigxifirst = sqrt((1-rhofirst^2)*sigomg^2);  
sigxisecond = sqrt((1-rhosecond^2)*sigomg^2);  

sigxilnw = sqrt((1-rholnw^2)*siglnw^2);  %/* SD of innovation in ln(wage) - set such that variance of ln(wage) is constant across time */

%/* Period 0 values of omega and ln(wage) */

omgdata0 = sigomg.*randn(n,1); 
lnwdata0 = siglnw.*randn(n,1);

%/* Period 1-b value of omega and Period 1 values of omega and ln(wage) */

omgdataminusb(:,1) = rhofirst*omgdata0 + sigxifirst.*randn(n,1);
omgdata(:,1) = rhosecond*omgdataminusb(:,1) + sigxisecond.*randn(n,1);
lnwdata(:,1) = rholnw*lnwdata0 + sigxilnw.*randn(n,1);

%/* Simulate values of omega and ln(wage) for rest of time periods */

i=1;
while i<=n
	s = 2;
 while s<=overallt
        omgdataminusb(i,s) = rhofirst*omgdata(i,s-1) + sigxifirst*randn(1,1);		
		omgdata(i,s) = rhosecond*omgdataminusb(i,s) + sigxisecond*randn(1,1);
		lnwdata(i,s) = rholnw*lnwdata(i,s-1) + sigxilnw*randn(1,1);

	s=s+1;
 end
	i=i+1;
end


lnkdata(:,1) = -100*ones(n,1);  %/* Initial ln(capital) level (very close to 0 capital) */ 

disc = 0.95;  %/* Discount rate for dynamic programming problem */
delta = 0.2;  %/* Depreciation rate for capital */

sigb = 0.6; %/* measures variation in capital adjustment cost across firms */
oneoverbiadj = exp(sigb.*randn(n,1)); %/* 1/capital adjustment cost (1/phi(i) in appendix) for each firm */

%/* Various components of long expression for optimal investment choice (end of ACF(2015) Appendix Section 7.3) */

squarebracketterm = (alphal^(alphal/(1-alphal)))*exp(0.5*alphal^2*sigoptl^2) - (alphal^(1/(1-alphal)))*exp(0.5*sigoptl^2); 
const1 = disc * (alphak/(1-alphal)) * (exp(alpha0)^(1/(1-alphal))) * squarebracketterm;
vec1 = (disc*(1-delta)).^seqa(0,1,100); 
vec2 = cumsumc(rholnw.^(2*seqa(0,1,100)));      
vec3 = sigxi^2 * [zeros(1,1);cumsumc(rho.^(2*seqa(0,1,99)))];  
expterm3 = exp( 0.5 * ((-alphal)/(1-alphal))^2 * sigxilnw^2 * vec2 );  
expterm4 = exp( 0.5 * (1/(1-alphal))^2 * rhosecond^2 * ((sigxifirst^2)*rho.^(2*seqa(0,1,100)) + vec3));
expterm5 = exp((1/(1-alphal))*(1/2)*sigxisecond^2);

%/* Compute optimal investment (and resulting capital stock) for all firms and all periods in simulation */

investmat = zeros(n,overallt);
i=1;
while i<=n

	s = 1;
 while s<=overallt;

		expterm1 = exp( (1/(1-alphal))*omgdata(i,s)*(rho.^seqa(1,1,100)));   %/* first term in exponent in second line */
		expterm2 = exp( ((-alphal)/(1-alphal))*lnwdata(i,s)*(rholnw.^seqa(1,1,100)) );   %/* second term in exponent in second line */

		investmat(i,s) = oneoverbiadj(i)*const1*expterm5*sum(vec1.*expterm1.*expterm2.*expterm3.*expterm4);
    if s>=2
       lnkdata(i,s) = log((1-delta)*exp(lnkdata(i,s-1)) + (1-0*rand(1,1))*investmat(i,s-1));
    end

		s=s+1;
 end
	i=i+1;
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%/* Now generate levels of labor input for all firms and time periods - note: this choice depends on omega(t-b) since labor chosen at t-b*/

i=1;
while i<=n
	s = 1;
  while s<=overallt
        lnldata(i,s) = ((sigxisecond^2)/2 + log(alphal) + alpha0 + rhosecond*omgdataminusb(i,s) - lnwdata(i,s) + lnpdata(i,s) + (alphak)*lnkdata(i,s))/(1-alphal);
		s=s+1;
  end
	i=i+1;
end

%/* Now add potential optimization error to ln(labor) - note: sigoptl defined at top of program */
truelnldata = lnldata;
lnldata = lnldata + sigoptl.*randn(n,overallt);

%/* Generate levels of Output and Materials */
lnydata = alpha0 + alphal*lnldata + alphak*lnkdata + omgdata + epsdata;
lnmdata = alpha0 + alphal*truelnldata + alphak*lnkdata + omgdata;  

DGP1 = [lnkdata(:,(starttime+1):(starttime+t)) lnldata(:,(starttime+1):(starttime+t)) lnmdata(:,(starttime+1):(starttime+t)) lnydata(:,(starttime+1):(starttime+t))]; 
cell2csv(DGP1, ['dgp',num2str(r)]); %export data to csv file

r = r+1;
end %end of simulation; replicate nrep times




