
%
% This program uses data from single_explore.do to
% estimate probits on plants closing down
%
%

clear variables

% Need to set path for Lesage's routines, etc.
 p = genpath('e:\copeland\Matlab_routines');
 addpath(p)

% Loading in data, creating variables
% 1.pin_sal_US 2.mrkt_price 3.prod 4.lag_levelprc 5.lag_levelsal 6.lag_levelinv
% 7.lag_levelinvsq 8.myindx 9.myindx_sq 10.prod_flg 11.inv_infer 12.cmpy 13.seg
% 14.myear 15.plt 16.mdl 17.indx 18.p_flg 19.days_supply 20.day_flg 21.I_flg

diary res_probitjul08c.txt; disp('diary on');

  % This is plant level data (as opposed to model level data)
  % Note that certain variables (e.g. model) are meaningless, see
  % single_explore.do to see construction of variables
  load data_moments_plt.csv
  data_moments=data_moments_plt;
  clear data_moments_plt
  
  sales  = data_moments(:,1); 
  price = data_moments(:,2); 
  prod = data_moments(:,3);
  lag_prc = data_moments(:,4); 
  lag_sal = data_moments(:,5); 
  lag_inv = data_moments(:,6); 
  lag_invsq = data_moments(:,7); 
  myindx = data_moments(:,8);
  myindx_sq = data_moments(:,9); 
  prod_flg = data_moments(:,10); 
  inv_infer = data_moments(:,11); 
  company = data_moments(:,12); 
  mrktseg = data_moments(:,13);
  myear = data_moments(:,14); 
  plant = data_moments(:,15); 
  model = data_moments(:,16); 
  month_indx = data_moments(:,17); 
  p_flg = data_moments(:,18);
  days_supply = data_moments(:,19); 
  day_flg = data_moments(:,20);
  I_flg = data_moments(:,21);
 

  clear data_moments
  disp('Data loaded');
  
%
% Creating fixed effects for plants
%
  nmbr_plt = max(plant);
  mat_plt = zeros(size(plant,1),nmbr_plt-1);
  for i=1:nmbr_plt-1
    indx = find(plant==i);
    mat_plt(indx,i)=1;
  end

%
% model year index squared
%
myindx2 = myindx.^2;


%
% Probit (Using Lesage's code)
%

  %detrending price, sales, and inventory - want differences from trend.
  indx = find(lag_prc>0 & lag_sal>0 & lag_inv>0 & sales>0 & price>0 & prod_flg>3 & p_flg<1);

  plt_dummy = mat_plt(indx,:);
 
  Y = price(indx);
  X = [ones(length(indx),1) myindx(indx) myindx2(indx)];
  p_res = ols(Y,X);
  
  Y = sales(indx);
  s_res = ols(Y,X);
 
  Y = inv_infer(indx);
  i_res = ols(Y,X);

  % creating lagged variables
  
%%  t_mat = [plant(indx) myear(indx) month_indx(indx) price(indx) sales(indx) inv_infer(indx) I_flg(indx) myindx(indx) myindx2(indx)];
  t_mat = [plant(indx) myear(indx) month_indx(indx) p_res.resid s_res.resid i_res.resid I_flg(indx) myindx(indx) myindx2(indx)];
  [m,n] = size(t_mat);
  sort_mat = sortrows(t_mat,[1 2 3]);
  X_lag1 = zeros(m,5+nmbr_plt-1);
  X_lag2 = zeros(m,8+nmbr_plt-1);
  Y_lag1 = zeros(m,1);
  Y_lag2 = zeros(m,1);
 
  j=1;
  for i=2:length(t_mat)
    if(t_mat(i,1)==t_mat(i-1,1) & t_mat(i,2)==t_mat(i-1,2))
      X_lag1(j,:)=[1 t_mat(i-1,4:6) t_mat(i,8) plt_dummy(i,:)];
      Y_lag1(j,1) = t_mat(i,7);
      j=j+1;
    end
  end 
  k=1;
  for i=3:length(t_mat)
    if(t_mat(i,1)==t_mat(i-2,1) & t_mat(i,2)==t_mat(i-2,2))
      X_lag2(k,:)=[1 t_mat(i-1,4:6) t_mat(i-2,4:6) t_mat(i,8) plt_dummy(i,:)];
      Y_lag2(k,1) = t_mat(i,7);
      k=k+1;
    end
  end

  

  X_lag1 = X_lag1(1:j-1,:);
  Y_lag1 = Y_lag1(1:j-1,:);
  X_lag2 = X_lag2(1:k-1,:);
  Y_lag2 = Y_lag2(1:k-1,:);
 
 
  % setting up the probit
  prb_res = probit(Y_lag1,X_lag1);
  prb_res2 = probit(Y_lag2,X_lag2);


  tmp1 = length(Y_lag1);  disp('One lag model, number of obs'); disp(tmp1);
  tmp1 = length(Y_lag2);  disp('Two lag model, number of obs'); disp(tmp1);

  % Compute marginal effects
%  tmp1 = mean(X)*prb_res.beta;
%  tmp2 = 1/sqrt(2*pi) * exp(-tmp1*tmp1/2);   %evaluate tmp1 w/ pdf standard normal
%  mrg_eff = prb_res.beta*tmp2;

  % Back out std err from tstats
  se = prb_res.beta./prb_res.tstat;

  disp('Probit results coef, se');
  disp([prb_res.beta(2:5) se(2:5)]);
  disp('Probit R-sqaured');
  disp(prb_res.rsqr);

  % Back out std err from tstats
  se = prb_res2.beta./prb_res2.tstat;

  disp('Probit results coef, se');
  disp([prb_res2.beta(2:8) se(2:8)]);
  disp('Probit R-sqaured');
  disp(prb_res2.rsqr);


  disp('Done');

diary off
