figure()
x07=[81	79	93	97	119	121	137	149	192];
x08=[141 94	126	119	119	150	149	168	233];
x09=[143 108 120 115 147 143 150 167 235];
x10=[99	104	95 119 137 143 145 169 229];
x11=[105 117 106 124 133 159 143 181 239];
x12=[107 107 105 143 127 152 176 177 256];
x13=[108 132 108 150 149 164 159 199 246];
x14=[115 113 133 133 132 153 164 164 233];
x15=[136 109 123 110 137 149 138 156 236];
cpi07=207.3;
cpi08=215.303;
cpi09=214.537;
cpi10=218.056;
cpi11=224.939;
cpi12=229.594;
cpi13=232.957;
cpi14=236.736;
cpi15=237.017;
nbins07=4;
nbins08=5;
nbins09=5;
% draw the first histogram
figure(1)
histogram(x07/cpi07,nbins07);
h_f1 = gca;    % get axis properties
d1 = h_f1;  % get Data properties under Histogram property set
% draw the second histogram
figure(2)
histogram(x08/cpi08,nbins08);
h_f2 = gca;
d2 = h_f2;
%draw the third histogram
figure(3)
histogram(x09/cpi09,nbins09);
h_f3 = gca;
d3 = h_f3;
% Drawing all histograms together
figure(4)
hist07=histogram(x07/cpi07,nbins07,'Displaystyle','stairs','Normalization','probability');
a07=hist07.Values;
b07=(hist07.BinEdges(1:max(size(hist07.BinEdges))-1)+hist07.BinEdges(2:max(size(hist07.BinEdges))))/2;
hold on;
hist08=histogram(x08/cpi08,nbins08,'Displaystyle','stairs','Normalization','probability');
a08=hist08.Values;
b08=(hist08.BinEdges(1:max(size(hist08.BinEdges))-1)+hist08.BinEdges(2:max(size(hist08.BinEdges))))/2;
hold on;
hist09=histogram(x09/cpi09,nbins09,'Displaystyle','stairs','Normalization','probability');
a09=hist09.Values;
b09=(hist09.BinEdges(1:max(size(hist09.BinEdges))-1)+hist09.BinEdges(2:max(size(hist09.BinEdges))))/2;
hold on;
hist10=histogram(x10/cpi10,nbins09,'Displaystyle','stairs','Normalization','probability');
a10=hist10.Values;
b10=(hist10.BinEdges(1:max(size(hist10.BinEdges))-1)+hist10.BinEdges(2:max(size(hist10.BinEdges))))/2;
hold on;
hist11=histogram(x11/cpi11,nbins09,'Displaystyle','stairs','Normalization','probability');
a11=hist11.Values;
b11=(hist11.BinEdges(1:max(size(hist11.BinEdges))-1)+hist11.BinEdges(2:max(size(hist11.BinEdges))))/2;
hold on;
hist12=histogram(x12/cpi12,nbins09,'Displaystyle','stairs','Normalization','probability');
a12=hist12.Values;
b12=(hist12.BinEdges(1:max(size(hist12.BinEdges))-1)+hist12.BinEdges(2:max(size(hist12.BinEdges))))/2;
hold on;;
hist13=histogram(x13/cpi13,nbins09,'Displaystyle','stairs','Normalization','probability');
a13=hist13.Values;
b13=(hist13.BinEdges(1:max(size(hist13.BinEdges))-1)+hist13.BinEdges(2:max(size(hist13.BinEdges))))/2;
hold on;
hist14=histogram(x14/cpi14,nbins09,'Displaystyle','stairs','Normalization','probability');
a14=hist14.Values;
b14=(hist14.BinEdges(1:max(size(hist14.BinEdges))-1)+hist14.BinEdges(2:max(size(hist14.BinEdges))))/2;
hold on;
hist15=histogram(x15/cpi15,nbins09,'Displaystyle','stairs','Normalization','probability');
a15=hist15.Values;
b15=(hist15.BinEdges(1:max(size(hist15.BinEdges))-1)+hist15.BinEdges(2:max(size(hist15.BinEdges))))/2;
%
title('US: Total Food Statistics');
legend('2007','2008','2009','2010','2011','2012','2013','2013','2014','2015');
xlabel('Income on Total Food'); ylabel('PDF of Income')
set(gca,'XLim')




figure(5)

xx = b07(1)-0.02:.01:b07(max(size(b07)))+0.02;
yy = spline(b07,a07,xx);
plot(b07,a07,'bo',xx,yy,'b-')
hold on

xx = b08(1)-0.02:.01:b08(max(size(b08)))+0.02;
yy = spline(b08,a08,xx);
plot(b08,a08,'ro',xx,yy,'r-')
hold on

xx = b09(1)-0.02:.01:b09(max(size(b09)))+0.02;
yy = spline(b09,a09,xx);
plot(b09,a09,'go',xx,yy,'g-')
hold on

xx = b10(1)-0.02:.01:b10(max(size(b10)))+0.02;
yy = spline(b10,a10,xx);
plot(b10,a10,'co',xx,yy,'c-')
hold on

xx = b11(1)-0.02:.01:b11(max(size(b11)))+0.02;
yy = spline(b11,a11,xx);
plot(b11,a11,'mo',xx,yy,'m-')
hold on

xx = b12(1)-0.02:.01:b12(max(size(b12)))+0.02;
yy = spline(b12,a12,xx);
plot(b12,a12,'bs',xx,yy,'b--')
hold on

xx = b13(1)-0.02:.01:b13(max(size(b13)))+0.02;
yy = spline(b13,a13,xx);
plot(b13,a13,'rs',xx,yy,'r--')
hold on

xx = b14(1)-0.02:.01:b14(max(size(b14)))+0.02;
yy = spline(b14,a14,xx);
plot(b14,a14,'gs',xx,yy,'g--')
hold on

xx = b15(1)-0.02:.01:b15(max(size(b15)))+0.02;
yy = spline(b15,a15,xx);
plot(b15,a15,'cs',xx,yy,'c--')
hold on


title('US: Total Food Statistics');
legend('2007','2007','2008','2008','2009','2009','2010','2010','2011','2011','2012','2012','2013','2013','2014','2014','2015','2015');
xlabel('Income on Total Food'); ylabel('PDF of Income')
set(gca,'XLim')
print(['us_histogram.png'],'-dpng')
print(['us_histogram.eps'],'-depsc')

figure(6)
plot(log(b07),log(a07),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1; 
a07(:,1)=[];
b07(:,1)=[];
x1=polyfit(log(b07),log(a07),1);
f1 = polyval(x1,log(b07));
plot(log(b07),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2007');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2007.png'],'-dpng')
print(['loglog_2007.eps'],'-depsc')

figure(7)
plot(log(b08),log(a08),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1;
a08(:,4)=[];
b08(:,4)=[];
a08(:,1)=[];
b08(:,1)=[];
x1=polyfit(log(b08),log(a08),1);
f1 = polyval(x1,log(b08));
plot(log(b08),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2008');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2008.png'],'-dpng')
print(['loglog_2008.eps'],'-depsc')

figure(8)
plot(log(b09),log(a09),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1; 
a09(:,4)=[];
b09(:,4)=[];
a09(:,1)=[];
b09(:,1)=[];
x1=polyfit(log(b09),log(a09),1);
f1 = polyval(x1,log(b09));
plot(log(b09),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2009');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2009.png'],'-dpng')
print(['loglog_2009.eps'],'-depsc')

figure(9)
plot(log(b10),log(a10),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1; 
a10(:,4)=[];
b10(:,4)=[];
a10(:,1)=[];
b10(:,1)=[];
x1=polyfit(log(b10),log(a10),1);
f1 = polyval(x1,log(b10));
plot(log(b10),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2010');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2010.png'],'-dpng')
print(['loglog_2010.eps'],'-depsc')

figure(10)
plot(log(b11),log(a11),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1; 
a11(:,4)=[];
b11(:,4)=[];
a11(:,1)=[];
b11(:,1)=[];
x1=polyfit(log(b11),log(a11),1);
f1 = polyval(x1,log(b11));
plot(log(b11),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2011');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2011.png'],'-dpng')
print(['loglog_2011.eps'],'-depsc')

figure(11)
plot(log(b12),log(a12),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1;
a12(:,4)=[];
b12(:,4)=[];
a12(:,1)=[];
b12(:,1)=[];
x1=polyfit(log(b12),log(a12),1);
f1 = polyval(x1,log(b12));
plot(log(b12),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2012');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2012.png'],'-dpng')
print(['loglog_2012.eps'],'-depsc')

figure(12)
hold on
plot(log(b13),log(a13),'bo')
h_f1 = gca;    % get axis properties
d1 = h_f1; 
a13(:,1)=[];
b13(:,1)=[];
x1=polyfit(log(b13),log(a13),1);
f1 = polyval(x1,log(b13));
plot(log(b13),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2013');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2013.png'],'-dpng')
print(['loglog_2013.eps'],'-depsc')

figure(13)
plot(log(b14),log(a14),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1; 
a14(:,4)=[];
b14(:,4)=[];
a14(:,1)=[];
b14(:,1)=[];
x1=polyfit(log(b14),log(a14),1);
f1 = polyval(x1,log(b14));
plot(log(b14),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2014');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2014.png'],'-dpng')
print(['loglog_2014.eps'],'-depsc')

figure(14)
plot(log(b15),log(a15),'bo')
hold on
h_f1 = gca;    % get axis properties
d1 = h_f1; 
a15(:,4)=[];
b15(:,4)=[];
a15(:,1)=[];
b15(:,1)=[];
x1=polyfit(log(b15),log(a15),1);
f1 = polyval(x1,log(b15));
plot(log(b15),f1)
legend(['slope= ',num2str(x1(1))]);
title('US PDF: Scaling Results');
legend('2015');
xlabel('Log(Income)'); ylabel('PDF of Income')
set(gca,'XLim')
print(['loglog_2015.png'],'-dpng')
print(['loglog_2015.eps'],'-depsc')