
% Time Series Disaggregation (TSD)

% clear all
% clc
% load EX3.mat

% n=number of HF data
% m=number LF data
% r=number of independant variables
% k=number of HF data between two consequent LF data
% beta[r*1]=parameters of regression
% Y[m*1]=LF dependant data
% X[n*1]=HF dependant data
% Z[n*r]=HF independant data

clear all
clc
load EX3

m=length(Y(:,1));
c=length(Y(1,:));
n=length(Z(:,1));
r=length(Z(1,:));
k=n/m;

% Constructing difference matrix D[n*n]

for i=1:c
    E(i,1)=1;
end

D=zeros(n,n);
D(1,1)=1;
for i=2:n
    for j=2:n
        if i==j
            D(i,j-1)=-1;
            D(i,j)=1;
        else
        end
    end
end

% Calculating A[n*n]

A=D'*D;

% Constructing aggregation matrix B[n*m]

B=zeros(n,m);
c0=0;
for j=1:m
    c0=c0+1;
    for i0=(c0-1)*k+1:c0*k
        B(i0,j)=1;
    end
end

% Calculating beta and X

PHI=(Z-inv(A)*B*inv(B'*inv(A)*B)*B'*Z)*inv(Z'*B*inv(B'*inv(A)*B)*B'*Z)*Z'*B*inv(B'*inv(A)*B)*Y+inv(A)*B*inv(B'*inv(A)*B)*Y;
SAI=(Z-inv(A)*B*inv(B'*inv(A)*B)*B'*Z)*inv(Z'*B*inv(B'*inv(A)*B)*B'*Z)*(Z'*B*inv(B'*inv(A)*B)*B'*inv(A)-Z')+(inv(A)*B*inv(B'*inv(A)*B)*B'*inv(A)-inv(A));
gama=pinv(SAI)*(C-PHI*E)*inv(E'*E);
beta=inv(Z'*B*inv(B'*inv(A)*B)*B'*Z)*Z'*B*inv(B'*inv(A)*B)*Y+inv(Z'*B*inv(B'*inv(A)*B)*B'*Z)*(Z'*B*inv(B'*inv(A)*B)*B'*inv(A)-Z')*gama*E';
lambda=inv(B'*inv(A)*B)*(B'*Z*beta-Y-B'*inv(A)*gama*E');
X1=PHI+SAI*gama*E';
X=PHI+(C-PHI*E)*inv(E'*E)*E';

sumX=0;
for q=1:c
    sumX=sumX+X(:,q);
end
DELTA=(C-sumX)'*(C-sumX);
RMSE1=sqrt(DELTA/n);

for q=1:c
    DELTA=DELTA+(Y(:,q)-B'*X(:,q))'*(Y(:,q)-B'*X(:,q));
end
RMSE2=sqrt(DELTA/(m*c));

% Correction of Y and C to satisfy YE-B'C=0.

% 1. Correction of Y.

Yc=B'*(PHI+(C-PHI*E)*inv(E'*E)*E');
Ylost1=(Y-Yc)*E;

% 2. Correction of C.

PHI1=(Z-inv(A)*B*inv(B'*inv(A)*B)*B'*Z)*inv(Z'*B*inv(B'*inv(A)*B)*B'*Z)*Z'*B*inv(B'*inv(A)*B)*B'*C+inv(A)*B*inv(B'*inv(A)*B)*B'*C;
SAI1=(Z-inv(A)*B*inv(B'*inv(A)*B)*B'*Z)*(inv(Z'*B*inv(B'*inv(A)*B)*B'*Z)*Z'*B*inv(B'*inv(A)*B)*B'*inv(A)-inv(Z'*B*inv(B'*inv(A)*B)*B'*Z)*Z')+inv(A)*B*inv(B'*inv(A)*B)*B'*inv(A)-inv(A);
Cc=PHI1+SAI1*pinv(B*B'*SAI1)*(B*Y*E-B*B'*PHI1);
Clost=Cc-C;
Ylost2=B'*Clost;
