%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Homework #6  Key
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Program  Name   : HW6_s08.m
%Chapter 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Question #9
%Answer:
%   Part a) The straight line model is y = Beta_0 + Beta_1 x + e. Base on 
%           the theory, we would expect the metal level to decrease as the 
%           distance increases. Thus, the slope should be negative.
%   Part b) Yes. As the distance from the plant increases, the 
%           concentration of clcium tends to decrease.
%   Part c) No. As the distance from the plant increases, the 
%           concentration of arsenic tends to increase.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question #12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nth Y X] = textread('blackbream.dat', '%s %f %f');
Xmean = mean(X);
Ymean = mean(Y);
ssize = size(X,1);
SS_xy = sum((X-Xmean) .* (Y-Ymean) );
SS_xx = sum((X-Xmean).^2);
beta_1 = SS_xy / SS_xx;                  %slope
beta_0 = Ymean - beta_1 * Xmean;         %y-intercept
%sprintf(' The regression line is y = %2.2f    %2.2f x ', beta_0 , beta_1  ) 
subplot(2,3,1) 
 plot(X, Y, 'x')
hold on
x_fit = [min(X), max(X)];
y_fit = beta_0 + beta_1 .* x_fit;
plot(x_fit, y_fit, '-');
   ylabel('Number of Strikes','FontName','Bold Arial', 'FontSize',10);
   xlabel('Age of Fish (days)','FontName','Bold Arial',...
   'FontSize',10);
   hold off;
   clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Question #12
%Answer:
%   Part a) y=Beta_0 + Beta_1 x + e
%   Part b) y=175.7033-0.8195x
%   Part c) Beta_0 does not have any pratical meaning since x=0 is not 
%                in the observed range.
%   Part d)   Beta_1 = - 0.8195 shows that for each addtional day of age, the mean number of 
%                  strikes is estimated  to decrease by 0.8195 strikes.   
%             
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
   
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question #13
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

[nth Y X] = textread('Ojuice.dat', '%s %f %f');
Xmean = mean(X);
Ymean = mean(Y);
ssize = size(X,1);
SS_xy = sum((X-Xmean) .* (Y-Ymean) );
SS_xx = sum((X-Xmean).^2);
beta_1 = SS_xy / SS_xx;                  %slope
beta_0 = Ymean - beta_1 * Xmean;         %y-intercept
%sprintf(' The regression line is y = %2.5f    %2.5f x ', beta_0 , beta_1  ) 
subplot(2,3,2) 
 plot(X, Y, 'x')
hold on
x_fit = [min(X), max(X)];
y_fit = beta_0 + beta_1 .* x_fit;
plot(x_fit, y_fit, '-');
   ylabel('Sweetness Index','FontName','Bold Arial', 'FontSize',10);
   xlabel('Pectin (ppm)','FontName','Bold Arial',...
   'FontSize',10);
 hold off;
 
 %%Part c
 y= beta_0 + beta_1.*300;
 clear;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Question #13
%Answer:
%   Part a) y=6.2521-.0023x
%   Part b) Since x=0 is not in the observed range, Beta_0= 6.25 has no interpretation other than being the 
%                 y-intercept. Beta_1 =-0.0023 shows that for each additional increase of 1 part per million of 
%                 pectin, the mean sweetness index is estimated to decreas by 0.0023.
%   Part c)   5.5589
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Question #20
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
[X Y] = textread('Heat.dat', ' %5.5f  %5.5f');
Xmean = mean(X);
Ymean = mean(Y);
ssize = size(X,1);
SS_xy = sum((X-Xmean) .* (Y-Ymean) );
SS_xx = sum((X-Xmean).^2);
beta_1 = SS_xy / SS_xx;                  %slope
beta_0 = Ymean - beta_1 * Xmean;         %y-intercept
%sprintf(' The regression line is y = %2.5f    %2.5f x ', beta_0 , beta_1  ) 
subplot(2,3,3) 
 plot(X, Y, 'x');
hold on
x_fit = [min(X), max(X)];
y_fit = beta_0 + beta_1 .* x_fit;
plot(x_fit, y_fit, '-');
   ylabel('Heat Transfer Enhancement','FontName','Bold Arial', 'FontSize',10);
   xlabel('Unflooded Area Ratio','FontName','Bold Arial',...
   'FontSize',10);
hold off;


 %%Part c
 SS_yy=sum((Y-Ymean).^2);
 SSE=(SS_yy - beta_1.*SS_xy);
 s=sqrt(SSE/(ssize-2));
 %s
 clear;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Question #20
%Answer:
%   Part a) y=0.2134+2.4264x
%   Part c) SSE=4.5311  s^2=0.2060
%   Part d) s=0.4538 .
%    We expect most of the observed y-values to lie within
% 2s (.9077) of their respective least squares predicted value since y is normally distributed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Question #26
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
[X Y] = textread('Heat.dat', ' %5.5f  %5.5f');
Xmean = mean(X);
Ymean = mean(Y);
ssize = size(X,1);
alpha=.05;
T=tinv(1-alpha/2,ssize-2);
SS_xx = sum((X-Xmean).^2);
SS_xy = sum((X-Xmean) .* (Y-Ymean) );
SS_xx = sum((X-Xmean).^2);
beta_1 = SS_xy / SS_xx;                  %slope
beta_0 = Ymean - beta_1 * Xmean;         %y-intercept
SS_yy=sum((Y-Ymean).^2);
 SSE=(SS_yy - beta_1.*SS_xy);
 s=sqrt(SSE/(ssize-2));
s_Beta=(s./(sqrt(SS_xx)));
lb=(beta_1-T.*s_Beta);
ub=(beta_1+T.*s_Beta);
%lb
%ub
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Question #26
%Answer: The 95% CI  lower_bound=1.9530  and the  upper_bound=2.8998
% Interpretation: We are 95% confident the mean heat transfer enhancement will increase from between
% 1.9530 and 2.8998 for each 1 unit increase in  unflooded area ratio.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Question #30
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




[Country Y X] = textread('GLOBRISK.dat', '%s %f %f');
Xmean = mean(X);
Ymean = mean(Y);
ssize = size(X,1);
SS_xy = sum((X-Xmean) .* (Y-Ymean) );
SS_xx = sum((X-Xmean).^2);
SS_yy = sum((Y-Ymean).^2);
beta_1 = SS_xy / SS_xx;                  %slope
beta_0 = Ymean - beta_1 * Xmean;         %y-intercept
%sprintf(' The regression line is y = %2.5f    %2.5f x ', beta_0 , beta_1  )

% plot
subplot(2, 3, 4)
plot(X, Y, 'x')
hold on
x_fit = [min(X), max(X)];
y_fit = beta_0 + beta_1 .* x_fit;
plot(x_fit, y_fit, '-');
xlabel('Credit Rating');
ylabel('Annualized Risk');
title('Original Regression (b)')
hold off;
SSE = SS_yy - beta_1 * SS_xy;
s = sqrt(SSE / (ssize - 2));
%alpha = 0.05;
df = ssize-2;
t_stat = beta_1 / (s / sqrt(SS_xx))
p_value=tcdf(-abs(t_stat), ssize-2) + (1- tcdf(abs(t_stat), ssize-2))
%t_cv = tinv(1-alpha/2, df)

clear;

% Use Microsolft Excel eliminate outliers (31.8, 87.0) and (32.6, 74.1) 
%Then, save as GREW.dat by choosing save as type "Text (Tab Delimited)"
% and regress again
%load globnew.dat;
[Country Y X] = textread('GRNEW.dat', '%s %f %f');
Xmean = mean(X);
Ymean = mean(Y);
ssize = size(X,1);
ssize
SS_xy = sum((X-Xmean) .* (Y-Ymean) )
SS_xx = sum((X-Xmean).^2)
SS_yy = sum((Y-Ymean).^2)
beta_1 = SS_xy / SS_xx                  %slope
beta_0 = Ymean - beta_1 * Xmean         %y-intercept
sprintf(' The regression line is y = %2.5f    %2.5f x ', beta_0 , beta_1  )
SSE = SS_yy - beta_1 * SS_xy
s = sqrt(SSE / (ssize - 2))
%alpha = 0.05;
%df = ssize-2;
t_stat = beta_1 / (s / sqrt(SS_xx))
%t_cv = tinv(1-alpha/2, df)
p_value=tcdf(-abs(t_stat), ssize-2) + (1- tcdf(abs(t_stat), ssize-2))


subplot(2,3,6);
scatter(X,Y,5,'x')
hold on;
int = min(X):1:max(X);
y_hat = beta_0 + (beta_1).*int;
plot(int,y_hat,'k');
   title('Outliers Eliminated','FontName','Bold Arial','FontSize',12);
   ylabel('Annulized Risk %','FontName','Bold Arial', 'FontSize',10);
   xlabel('Average Credit Raiting','FontName','Bold Arial',...
   'FontSize',10);
 hold off;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Question #30
%Answer:
%   Part a) Beta_0=57.755         Beta_1=-0.3996
%   Part b) y=57.7551-0.3996x
%   Part c) To determine if country credit risk contributes information for
%   the prediction of market volatility, we test H_0: beta_1=0 vs beta_1
%   ~=0. The observed test statistic t=-4.3665 and p-value=0.0000938. So
%   H_0 is rejected.
%   Part d) Delet outliers (31.8,87.0) and (32.6, 74.1). 
%   Part e) After deleting outliers, we repeat previous steps to get y=47.0323-0.2633x
%   reconstruct picture and test. The results do not change dramatically.  The observed 
%   t=-3.7583 and p-valued =
%   0.0006065. H_0 is rejected
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
 
 
