%gaseprice.m
%Gas price per 1000 ft3 1/81 - 12/13
load gas;
n = length(gas);
t = 1:n; t=t';
yrs = 1981: 1/12:2014 - 1/12;
figure(1)
plot(yrs,gas)
title('Monthly Residential Gas Price (per 1000 ft3) 1981-2013')
xlabel('Year')
ylabel('$')
grid 'minor'
pause
%=================================================
% Various Definitions of E[X(t)] = mu(t):
% DEFINITION 1:
% mu(t)defined as the slowly changing trend:
% mu(t)HAT = 12-month moving average:
gavg=zeros(n,1);
for k = 13:n
    gavg(k)=mean(gas(k-11:k));
end
gavg(1:12)=gavg(13)*ones(12,1);
%DEFINITION 2:
% mu(t) = mu(k) for the kth year:
g12avg = [];
for k = 1:fix(n/12)
    g12(k)=mean(gas(1+(k-1)*12:k*12));
    g12avg = [g12avg , g12(k)*ones(1,12)];
end
g12avg = g12avg';
dg = gas - gavg;
w=2*pi/12; % Annual frequency (rad/sec)
a1=2*cos(w); a2=-1;
dghat=dg;
for k=3:n
    dghat(k)= a1*dg(k-1) + a2*dg(k-2);
end
%------------------------------------
ddg = dg - dghat;
D = [gas,gavg,g12avg,dg,dg12,dghat,ddg];
figure(1)
plot(yrs,D)
ylabel('$')
grid 'minor'
legend('gas','MA','12moA','dg','dg12','dghat','ddg')
gpr= gavg + dghat; %LEVEL 1 PREDICTOR
pause
%=====================================================
figure(2)
plot(yrs,[gas,gpr])
ylabel('$')
grid 'minor'
legend('gas','1-month ahead LEVEL 0 prediction')
pause
no=25;
dgpr = gas(no:n) - gpr(no:n);
figure(3)
plot(yrs(no:n),dgpr)
ylabel('$')
grid 'minor'
legend('LEVEL 0 prediction residual')
msedgpr = mean(dgpr.^2)
pause
%======================================================
ffto=abs(fft(dgpr)); ffto = ffto(1:nn/2);
nn = length(dgpr);
fsamp = 1; %samples per month
df = fsamp/nn;
fvec = 0:df:fsamp-df; fvec=fvec(1:nn/2)';
figure(3)
plot(fvec,ffto)
title('|FT(dgpr)|')
xlabel('Frequency (cycles/month)')
grid 'minor'
pause
%=======================================================
[fftomax,I]=max(ffto);
fmax = df*I; % Frequency of max residual
Tmax = (1/fmax); %Period of max residual = 6 months
g6pr = zeros(n-no+1,1);
a11 = cos(2*pi/Tmax); a22 = -1;
for k = 3:n-no+1
    g6pr(k) = a11*dgpr(k-1) + a22*dgpr(k-2);
end
dgpr6 = dgpr - g6pr;
figure(4)
plot(yrs(no:n),[dgpr , dgpr6])
legend('pr','ddd')
title('Level 0 (dgpr) and LEVEL 1 (dgpr6) Prediction Residuals')
grid 'minor'
pause
figure(3)
hold on
dddft=abs(fft(dgpr6)); dddft=dddft(1:nn/2);
plot(fvec,dddft,'r')
title('|FT(dgpr)|(blue) and |FT(dgpr6)| (red)')
msedgpr6 = mean(dgpr6.^2)


