function diagfore=diagforecast(data,p,q)
k=size(data,2);
[parameters, Ht]=diagonal_bekk_mvgarch(data,p,q, []);
k2=k*(k+1)/2;
C=parameters(1:(k2));
A=parameters(k2+1:k2+k*p);
B=parameters(k2+k*p+1:k2+k*p+k*q);
tempA=zeros(k,k,p);
tempB=zeros(k,k,p);
for i=1:p
    tempA(:,:,i)=diag(A((k*(i-1)+1):(k*i)));
end
for i=1:q
    tempB(:,:,i)=diag(B((k*(i-1)+1):(k*i)));
end
A=tempA;
B=tempB;

C=ivech(C);
C=tril(C);

Htt=Ht(:,:,size(Ht,3));
epsi=data(size(data,1),:);
epsilon=epsi'*epsi;
diagfore=C*C'+B*Htt*B'+A*epsilon*A';