%Simulate Data and Predict
%y=a + bx + cx^2
n=20;
a=5.0;
b=2.25;

%simulate
x=2.5*[1:n]';
yactual=[a+b*x];

sigma=5;
e=[normrnd(0,sigma,n,1)];
y=yactual+e;

%estimate
xmat=[[ones(1,n)]',x];
xdim=size(xmat);
p=xdim(2);

est=inv(xmat'*xmat)*(xmat'*y);
yfit=xmat*est;
sigma_est=(y-yfit)'*(y-yfit)/(n-p);
var_est=sigma_est*inv(xmat'*xmat);
se_est=sqrt(diag(var_est));
tstat=est./se_est;


%Prediction is yfit
%Prediction variance

xnew=[1:0.1:(2.5*n)]';

xnewmat=[[ones(1,length(xnew))]',xnew];
ynew=xnewmat*est;

hold on
plot(x,y,'Color','red','Marker','*','LineStyle','none');
plot(xnew,ynew,'Color','green');
axis([0 50 0 100]);

var_ynew=diag(xnewmat * (var_est * xnewmat'));

var_ynew_pred=diag(xnewmat * (var_est * xnewmat'))+sigma_est;

%Prediction CI
p_upper=ynew+1.96*sqrt(var_ynew);
p_lower=ynew-1.96*sqrt(var_ynew);

plot(xnew,p_upper,'--g',xnew,p_lower,'--g');

%Prediction CI
o_upper=ynew+1.96*sqrt(var_ynew_pred);
o_lower=ynew-1.96*sqrt(var_ynew_pred);

plot(xnew,o_upper,':g',xnew,o_lower,':g');

legend('Obs.','Pred.','Pred. CI','Obs. CI',2);
hold off;





