clear all
close all;
clc;
tic;

%global y z pi ai h W W1 i

load full.txt;
d=full;

load psiconf.txt;
psiconfi=psiconf;
psi5=psiconfi(:,4);
psi95=psiconfi(:,2);
psihati=psiconfi(:,7);
nrconf=psiconfi(:,5);
phi5=psiconfi(:,3);
phi95=psiconfi(:,1);
phihati=psiconfi(:,8);
naconf=psiconfi(:,6);

Pop2010=d(:,1);
under5	=d(:,2);
under18	=d(:,3);
over65	=d(:,4);
Female	=d(:,5);
Whitee	=d(:,6);
Black	=d(:,7);
AmeIndianAlaska	=d(:,8);
Asian	=d(:,9);
Two_or_More_Races	=d(:,10);
Hispanic_or_Latino	=d(:,11);
WhitenotH	=d(:,12);
same_house=d(:,13);
Foreignborn	=d(:,14);
otherEnglish =d(:,15);
Highs =d(:,16);
Bachelor =d(:,17);
Veterans =d(:,18);
travel_time	=d(:,19);
Housing_units =d(:,20);
Homeownership =d(:,21);
Housing_units_in_multi_unit	=d(:,22);
Median_value_of_owner_occupied_h	=d(:,23);
Households	=d(:,24);
per_household	=d(:,25);
Percapomoneyinc	=d(:,26);
Medianhhinc	=d(:,27);
below_poverty	=d(:,28);
Retail	=d(:,29);
Retailpercap	=d(:,30);
Land_area	=d(:,31);
per_sqmile	=d(:,32);
sincome	=d(:,33);
sprint	=d(:,34);
m_inches	=d(:,35);
m_dollars	=d(:,36);
m_unit	=d(:,37);
adrate	=d(:,38);
m_peru	=d(:,39);
weekday_circ	=d(:,40);
weekend_circ	=d(:,41);
price_yearly	=d(:,42);
income	=d(:,43);
wageinprint	=d(:,44);
avr_pageno	=d(:,45);
pop2011	=d(:,46);
wprinting	=d(:,47);
earnpaper	=d(:,48);
wpaper	=d(:,49);
gdp_growth	=d(:,50);
nr	=d(:,51);
na	=d(:,52);
loginc	=d(:,53);
logprint	=d(:,54);
logpop	=d(:,55);
extra	=d(:,56);
totalci	=d(:,57);
na2	=d(:,58);
days	=d(:,59);
dprice	=d(:,60);
logminc	=d(:,61);
logretail=log(Retail);
nr2=weekday_circ./Pop2010;
logpsq=log(per_sqmile);

b=500; %number of bootstrap replications

y=nr;
z1=na2;
z=na;

n=length(y);

D=npest_instd(y,z1,dprice,Land_area,Retail);
hy1=D(:,1);
phi1=D(:,2);

figure('Units','inches','Position',[0 0 9 4.5]);
figure(1);
subplot(1,2,1);
%set(gca,'FontSize', 12);
plot(y,hy1,'r+');
hold on;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^r$', 'interpreter', 'latex','FontSize',11);
subplot(1,2,2);
%set(gca,'FontSize', 12);
plot(z1,phi1,'r+');
hold on;
xlabel('Share of advertisers, $N^a$','interpreter', 'latex','FontSize',11); ylabel('Network effect on readers, $\varphi$','interpreter', 'latex','FontSize',11);
print -deps revhy.eps


DA=npest_instadv(z1,y,adrate,below_poverty,wageinprint);
hz1=DA(:,1);
psi1=DA(:,2);

figure('Units','inches','Position',[0 0 9 4.5]);
figure(2);
subplot(1,2,1);
plot(z1,hz1,'r+');
hold on;
xlabel('Share of advertisers $N^a$','interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^a$', 'interpreter', 'latex','FontSize',11);
subplot(1,2,2);
plot(y,psi1,'r+');
hold on;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Network effect on advertisers, $\psi$','interpreter', 'latex','FontSize',11);
print -deps revhz.eps

yy=sort(y);
hy11=sort(hy1,'descend');
zz=sort(z1);
hz11=sort(hz1,'descend');

figure('Units','inches','Position',[0 0 9 4.5]);
figure(3); %figure for monotonicity imposed by rearrangement
subplot(1,2,1);
plot(yy,hy11,'r+');
hold on;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^r$', 'interpreter', 'latex','FontSize',11);
subplot(1,2,2);
plot(z1,phi1,'r+');
hold on;
xlabel('Share of advertisers, $N^a$','interpreter', 'latex','FontSize',11); ylabel('Network effect on readers, $\varphi$','interpreter', 'latex','FontSize',11);
print -deps revhyr.eps

figure('Units','inches','Position',[0 0 9 4.5]);
figure(4); %figure for monotonicity imposed by rearrangement
subplot(1,2,1);
plot(zz,hz11,'r+');
hold on;
xlabel('Share of advertisers $N^a$','interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^a$', 'interpreter', 'latex','FontSize',11);
subplot(1,2,2);
plot(y,psi1,'r+');
hold on;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Network effect on advertisers, $\psi$','interpreter', 'latex','FontSize',11);
print -deps revhzr.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%robustness%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D1=npest_instd(y,z1,dprice,wageinprint,Retailpercap);
hyrob=D1(:,1);
phi1rob=D1(:,2);

DA1=npest_instadv(z1,y,adrate,Highs,wageinprint);
hzrob=DA1(:,1);
psi1rob=DA1(:,2);

hyrob2=sort(hyrob,'descend');
hzrob2=sort(hzrob,'descend');

figure('Units','inches','Position',[0 0 9 4.5]);
figure(5);
subplot(1,2,1);
plot(yy,hyrob2,'r+');
hold on;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^r$', 'interpreter', 'latex','FontSize',11);
subplot(1,2,2);
plot(z1,phi1rob,'r+');
hold on;
xlabel('Share of advertisers, $N^a$','interpreter', 'latex','FontSize',11); ylabel('Network effect on readers, $\varphi$','interpreter', 'latex','FontSize',11);
print -deps robustreader.eps

figure('Units','inches','Position',[0 0 9 4.5]);
figure(6);
subplot(1,2,1);
plot(zz,hzrob2,'r+');
hold on;
xlabel('Share of advertisers $N^a$','interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^a$', 'interpreter', 'latex','FontSize',11);
subplot(1,2,2);
plot(y,psi1rob,'r+');
hold on;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Network effect on advertisers, $\psi$','interpreter', 'latex','FontSize',11);
print -deps robustadv.eps



%%bootstrap for confidence intervals%%%%%%%%%%%%%%%%%%%%%%%
p=dprice;
pi=Land_area;
nai=Retail;
a=adrate;
ai=wageinprint;
ci=below_poverty;
dat=[y z1 p a pi ai nai ci];
[bts sam]=bootstrp(b,[],dat);
dn=zeros(n,8);
for k=1:b;
for i=1:n;
    dn(i,:)=dat(sam(i,k),:);
end;
yb=dn(:,1);
zb=dn(:,2);
pb=dn(:,3);
ab=dn(:,4);
pib=dn(:,5);
aib=dn(:,6);
naib=dn(:,7);
cib=dn(:,8);

DB=npest_instd(yb,zb,pb,pib,naib);
hy1b=DB(:,1);
phi1b=DB(:,2);

DBA=npest_instadv(zb,yb,ab,aib,cib);
hz1b=DBA(:,1);
psi1b=DBA(:,2);

yyb=sort(yb);
hyb11=sort(hy1b,'descend');
zzb=sort(zb);
hzb11=sort(hz1b,'descend');


for i=1:n;
    kk1=hy1b(yb==y(i));
    kk2=phi1b(yb==y(i));
    kk3=hz1b(yb==y(i));
    kk4=psi1b(yb==y(i));
    kk1r=hyb11(yyb==yy(i));
    kk3r=hzb11(yyb==yy(i));
    
if size(kk1)==[0,1]
    ahy(i,1)=0;
    aphi(i,1)=0;
    ahz(i,1)=0; 
    apsi(i,1)=0;
else
 ahy(i,1)=kk1(1,1);
 aphi(i,1)=kk2(1,1);
 ahz(i,1)=kk3(1,1);
 apsi(i,1)=kk4(1,1);
end;

 if size(kk1r)==[0,1]
     ahyr(i,1)=0;
     ahzr(i,1)=0;
 else
     ahyr(i,1)=kk1r(1,1);
     ahzr(i,1)=kk3r(1,1);
 end;
end;

kh1(:,k)=ahy;
kphi(:,k)=aphi;
kh2(:,k)=ahz;
kpsi(:,k)=apsi;
kh1r(:,k)=ahyr;
kh2r(:,k)=ahzr;

end;



r1=zeros(1,n);
r1r=zeros(1,n);
r2=zeros(1,n);
r3=zeros(1,n);
r4=zeros(1,n);
r3r=zeros(1,n);
for i=1:n;
    [l1,j1,v1]=find(kh1(i,:)); 
    r1(i)=std(v1);
    K1=sqrt(((v1-hy1(i)).^2)/(r1(i)^2));
    k1alphal(i)=quantile(K1,0.05);
    k1alphau(i)=quantile(K1,0.95);
    
    [l2,j2,v2]=find(kphi(i,:)); 
    r2(i)=std(v2);
    K2=sqrt(((v2-phi1(i)).^2)/(r2(i)^2));
    k2alphal(i)=quantile(K2,0.05);
    k2alphau(i)=quantile(K2,0.95);
    
    [l3,j3,v3]=find(kh2(i,:)); 
    r3(i)=std(v3);
    K3=sqrt(((v3-hz1(i)).^2)/(r3(i)^2));
    k3alphal(i)=quantile(K3,0.05);
    k3alphau(i)=quantile(K3,0.95);
    
    [l4,j4,v4]=find(kpsi(i,:)); 
    r4(i)=std(v4);
    K4=sqrt(((v4-psi1(i)).^2)/(r4(i)^2));
    k4alphal(i)=quantile(K4,0.05);
    k4alphau(i)=quantile(K4,0.95);
    
    [l1r,j1r,v1r]=find(kh1r(i,:)); 
    r1r(i)=std(v1r);
    K1r=sqrt(((v1r-hy11(i)).^2)/(r1r(i)^2));
    k1ralphal(i)=quantile(K1r,0.05);
    k1ralphau(i)=quantile(K1r,0.95);
    
    [l3r,j3r,v3r]=find(kh2r(i,:)); 
    r3r(i)=std(v3r);
    K3r=sqrt(((v3r-hz11(i)).^2)/(r3r(i)^2));
    k3ralphal(i)=quantile(K3r,0.05);
    k3ralphau(i)=quantile(K3r,0.95);
    
end;

ubhy=hy1'+k1alphau.*r1;
lbhy=hy1'-k1alphau.*r1;

ubphi=phi1'+k2alphau.*r2;
lbphi=phi1'-k2alphau.*r2;

ubhz=hz1'+k3alphau.*r3;
lbhz=hz1'-k3alphau.*r3;

ubpsi=psi1'+k4alphau.*r4;
lbpsi=psi1'-k4alphau.*r4;

ubrhy=hy11'+k1ralphau.*r1r;
lbrhy=hy11'-k1ralphau.*r1r;

ubrhz=hz11'+k3ralphau.*r3r;
lbrhz=hz11'-k3ralphau.*r3r;




figure('Units','inches','Position',[0 0 9 4.5]);
figure(7);
plot(z1,ubphi,'b+');
hold on;
plot(z1,lbphi,'b+');
hold on;
plot(z1,phi1,'r.');
hold off;
xlabel('Share of advertisers, $N^a$','interpreter', 'latex','FontSize',11); ylabel('Network effect on readers, $\varphi$','interpreter', 'latex','FontSize',11);
print -deps revphi.eps

figure('Units','inches','Position',[0 0 9 4.5]);
figure(8);
plot(y,ubpsi,'b+');
hold on;
plot(y,lbpsi,'b+');
hold on;
plot(y,psi1,'r.');
hold off;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Network effect on advertisers, $\psi$','interpreter', 'latex','FontSize',11);
print -deps revpsi.eps


figure('Units','inches','Position',[0 0 9 4.5]);
figure(9);
plot(yy(5:112),ubrhy(5:112),'b+');
hold on;
plot(yy(5:112),lbrhy(5:112),'b+');
hold on;
plot(yy(5:112),hy11(5:112),'r.');
hold off;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^r$', 'interpreter', 'latex','FontSize',11);
print -deps revhycir.eps

figure(10);
plot(y,ubhy,'b.');
hold on;
plot(y,lbhy,'b.');
hold on;
plot(y,hy1,'r.');
hold off;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^r$', 'interpreter', 'latex','FontSize',11);

figure('Units','inches','Position',[0 0 9 4.5]);
figure(11);
plot(z1,ubhz,'b.');
hold on;
plot(z1,lbhz,'b.');
hold on;
plot(z1,hz1,'r.');
hold off;
xlabel('Share of advertisers $N^a$','interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^a$', 'interpreter', 'latex','FontSize',11);



figure('Units','inches','Position',[0 0 9 4.5]);
figure(12);
plot(zz(3:112),ubrhz(3:112),'b+');
hold on;
plot(zz(3:112),lbrhz(3:112),'b+');
hold on;
plot(zz(3:112),hz11(3:112),'r.');
hold off;
xlabel('Share of advertisers $N^a$','interpreter', 'latex','FontSize',11); ylabel('Inverse of survival, $H^a$', 'interpreter', 'latex','FontSize',11);
print -deps revhzcir.eps

figure('Units','inches','Position',[0 0 9 4.5]);
figure(13);
subplot(1,2,1);
plot(y,z1,'.');
hold on;
xlabel('Share of readers,$N^r$','interpreter', 'latex','FontSize',11); ylabel('Share of advertisers,$N^a$','interpreter', 'latex','FontSize',11);
subplot(1,2,2);
plot(weekday_circ,m_inches,'.');
xlabel('Circulation','interpreter', 'latex','FontSize',11);ylabel('Advertising','interpreter', 'latex','FontSize',11);
hold off;
print -deps scatters.eps

phipara=10.5788-(48.3402*z1)+(96.7471*z1.^2)-(62.0625*z1.^3);
psipara=1.0588-(162.191*y)+(3912.491*y.^2);



figure('Units','inches','Position',[0 0 9 4.5]);
figure(14);
plot(z1,phipara,'b*');
hold on;
plot(z1,phi1,'r.');
hold on;
plot(z1,ubphi,'g+');
hold on;
plot(z1,lbphi,'g+');
hold off;
xlabel('Share of advertisers, $N^a$','interpreter', 'latex','FontSize',11); ylabel('Network effect on readers, $\varphi$','interpreter', 'latex','FontSize',11);
ll=legend('$\hat{\varphi}$ parametric', '$\hat{\varphi}$ nonparametric','nonparametric uCI','nonparametric lCI')
 set(ll,'Interpreter','latex','FontSize',12);
print -deps paranonparaphi.eps



figure('Units','inches','Position',[0 0 9 4.5]);
figure(15);
plot(y,psipara,'b*');
hold on;
plot(y,psi1,'r.');
hold on;
plot(y,ubpsi,'g+');
hold on;
plot(y,lbpsi,'g+');
hold off;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Network effect on advertisers, $\psi$','interpreter', 'latex','FontSize',11);
ll2=legend('$\hat{\psi}$ parametric', '$\hat{\psi}$ nonparametric','nonparametric uCI','nonparametric lCI')
set(ll2,'Interpreter','latex','FontSize',12);
print -deps paranonpara.eps

figure('Units','inches','Position',[0 0 9 4.5]);
figure(16);
plot(nrconf,psi5,'b+');
hold on;
plot(nrconf,psi95,'b+');
hold on;
plot(nrconf,psihati,'r.');
hold off;
xlabel('Share of readers, $N^r$', 'interpreter', 'latex','FontSize',11); ylabel('Network effect on advertisers, $\psi$','interpreter', 'latex','FontSize',11);
print -deps psipara.eps

figure('Units','inches','Position',[0 0 9 4.5]);
figure(17);
plot(naconf,phi5,'b+');
hold on;
plot(naconf,phi95,'b+');
hold on;
plot(naconf,phihati,'r.');
hold off;
xlabel('Share of advertisers, $N^a$','interpreter', 'latex','FontSize',11); ylabel('Network effect on readers, $\varphi$','interpreter', 'latex','FontSize',11);
print -deps phipara.eps

duree=toc ;
disp(['Duration = ',num2str(duree),' secondes']);