% Calculate two-first order coefficents for different parameterization of the
% storage model Perform distribution of first- and second-order autocorrelation
% coefficents on all consecutive subsamples of length 100

% In this model, shocks are generated by a normal with mean and standard
% deviation set to 1 and 0.05. The unconditional mean price is equal to 1.
% Storage cost and demand elasticity are varied to analyze their effect
% on serial correlation. We use the Proposition 1 of Deaton & Laroque (1996)
% to rescale the distritbuion to a truncated unit normal. This leads to the 
% following constraints on the parameters:
% a = 1
% Elasticity = 0.05/b
% Or with a truncated normal with mean 1 and standard deviation 0.05:
% a+b = 1
% Elasticity = 1/b

Elasticity0     = -0.05;
StorageCost0    = 0;
a = 1;

[~,~,~,~,~,~,~,options,~,N,mina,maxa,r,delta] = OverallInitialization;

%% Define the grid of interpolation (equally space points)
interp.s           = linspace(mina,maxa,N)';
interp.PriceInterp = [];

%% Define the shocks structure by a discretization of the gaussian density
truncated_normal_rule(3,10,0,1,-5,5,fullfile('..','Results','GaussianQuadrature'))
model.shocks.e  = load(fullfile('..','Results','GaussianQuadrature_x.txt'));
model.shocks.w  = load(fullfile('..','Results','GaussianQuadrature_w.txt'));
TruncatedNormal = truncate(makedist('Normal'),-5,5);
model.funrand   = @(nrep,nper) TruncatedNormal.random(nrep,nper);

%% Number of periods and replications
Nsample         = 2E4;                % Number of replications in the MC experiment
burninper       = 30;
samplelen       = 100;
nper            = 49+samplelen+burninper;
A0              = ones(Nsample,1); % Initial vector of availability

disp('Number of samples of length T on the asymptotic distribution')
disp(Nsample*(nper+1-burninper-samplelen))

%% Storage costs
StoCost         = linspace(0,0.2,100)';
corstatk        = NaN(Nsample*(nper+1-burninper-samplelen),2,length(StoCost));
Qcorstatk       = NaN(length(StoCost),2,3);

b = 0.05/Elasticity0;

for i = 1:length(StoCost)
  k             = StoCost(i);
  model.params  = [a b delta k r];
  interp        = SolveStorageDL(model,interp,options);
  rng(options.seed)
  [~,Xsim]      = SimulStorage(model,interp,A0,nper);
  P             = squeeze(Xsim(:,2,(burninper+1):end))';

  for t=1:(nper+1-burninper-samplelen)
    corstatk((1+(t-1)*Nsample):(t*Nsample),:,i) = autocor(P(t:(t+samplelen-1),:));
  end
  Qcorstatk(i,:,:)  = quantile(corstatk(:,:,i),[0.05 0.5 0.95])';
end

%% Demand elasticity
eta       = linspace(-0.4,-0.005,300)';
b         = 0.05./eta;
corstatb  = NaN(Nsample*(nper+1-burninper-samplelen),2,length(b));
Qcorstatb = NaN(length(b),2,3);

for i = 1:length(b)
  model.params  = [a b(i) delta 0 r];
  interp        = SolveStorageDL(model,interp,options);
  rng(options.seed)
  [~,Xsim]      = SimulStorage(model,interp,A0,nper);
  P             = squeeze(Xsim(:,2,(burninper+1):end))';

  for t=1:(nper+1-burninper-samplelen)
    corstatb((1+(t-1)*Nsample):(t*Nsample),:,i) = autocor(P(t:(t+samplelen-1),:));
  end
  Qcorstatb(i,:,:)  = quantile(corstatb(:,:,i),[0.05 0.5 0.95])';
end

%% Save results
save(fullfile('..','Results','SimulatedAutocorr.mat'),...
     'StoCost','eta','b','Qcorstatk','Qcorstatb')
