%PROGRAM NAME: ex9_13.m   (9/13/19)
%Create pressure blowdown data
%Create mean pressure mup=b1+b2*t: 
t=0:.01:20-.01;
b1=5; b2=-1/20;
mup=b1+b2*t;
%Create dp data:
nt=length(t);
a1=0.9; a2=0.2; %Range of AR(1) BW parameters
da=(a1-a2)/(nt-2);
a=a1:-da:a2-da; %AR(1) BW parameter vector
std1=0.5; std2=0.1; %Range of std_dp
dstd=(std1-std2)/(nt-2);
std_dp=std1:-dstd:std2-dstd;
dp=zeros(1,nt);
dp(1)=normrnd(0,std_dp(1));
v=zeros(1,nt);
for k=2:nt
    v(k)=normrnd(0,std_dp(k)*sqrt(1-a(k)^2));
    dp(k)=a(k)*dp(k-1)+v(k);
end
p=mup+dp;
figure(31)
plot(t,p)
title('Turbine Pressure Blow-Down Data')
xlabel('time (sec)')
ylabel('pressure (kpsi)')
grid
savefile='pressuredata.mat';
save(savefile,'t','p')
%=========================
load pressuredata
%This will result in the arrays 't' and 'p' being in your workspace
%Mean MODEL: phat=b1+b2*T:
Cov=cov(t',p');
b2hat=Cov(1,2)/Cov(1,1);
b1hat=mean(p)-b2hat*mean(t);
muphat=b1hat+b2hat*t;
hold on
plot(t,muphat,'r')
[b1 b2;b1hat b2hat]
%Residual Analysis:
Dp=p-muphat;
figure(32)
plot(t,Dp)
title('Residual Pressure Data')
xlabel('time (sec)')
ylabel('pressure (kpsi)')
grid
%Use a T-point MA to estimate std_Dp:
T=101;
std_Dp=zeros(1,nt);
for k=1+(T-1)/2 : nt-1-(T-1)/2
    std_Dp(k)=std(Dp(k-(T-1)/2 : k+(T-1)/2));
end
std_Dp(1:(T-1)/2)=std_Dp(1+(T-1)/2);
std_Dp(nt-(T-1)/2:nt)=std_Dp(nt-1-(T-1)/2);
figure(33)
plot(t,[std_dp;std_Dp])
title('Plots of True and Estimated std_P')
xlabel('time (sec)')
grid
%Scale residual pressur by std_Dp:
SDp=dp./std_Dp;
figure(34)
plot(t,SDp)
title('Scaled Residual Pressure Data')
xlabel('time (sec)')
ylabel('pressure (kpsi)')
grid
%Use a T-point Sliding AR(1) Model to estimate BW parameter a(t):

ahat=zeros(1,nt);
for k=1+(T-1)/2 : nt-1-(T-1)/2
    R=xcorr((SDp(k-(T-1)/2 : k+(T-1)/2)),1);
    ahat(k)=R(1)/R(2);
end
ahat(1:(T-1)/2)=ahat(1+(T-1)/2);
ahat(nt-(T-1)/2:nt)=ahat(nt-1-(T-1)/2);
figure(35)
plot(t,[a;ahat])
title('Plots of True and Estimated AR(1) BW Parameter a(t)')
xlabel('time (sec)')
%Test for whiteness of the scaled residual prediction error:
SDphat=ahat(2:nt).*SDp(1:nt-1);
SDperr=SDp(2:nt)-SDphat;
figure(36)
plot(t(2:nt),SDperr)
title('Plot of Scaled Residual Prediction Error')
xlabel('time (sec)')
figure(37)
errcorr=xcorr(SDperr,10);
tau=-10:10;
plot(tau,errcorr)
title(' Scaled Residual Pressure Prediction Error Sample Autocorrelation')
xlabel(' correlation lag')
grid

%Construct estimate of the original data: