Code covered by the BSD License

# Risk and Asset Allocation

### Attilio Meucci (view profile)

16 Nov 2005 (Updated )

Software for quantitative portfolio and risk management

[E_EM, S_EM, Y, CountLoop]=EM(X)
```function [E_EM, S_EM, Y, CountLoop]=EM(X)
% this function implements the Expectation-Maximization (EM) algorithm to recover
% missing observations in a time series
% see "Risk and Asset Allocation"- Springer (2005), by A. Meucci

[T,N]=size(X);

% E-M initialization
X_Init=[];
for t=1:T
if isempty(X(isnan(X(t,:))))
X_Init=[X_Init
X(t,:)];
end
end
M=mean(X_Init)';
S=cov(X_Init);

Tolerance=10^(-6)*mean(  [M
sqrt(diag(S))]  );

% E-M loop
Convergence=0;
CountLoop=0;
Y=X;
while ~Convergence
CountLoop=CountLoop+1;

% Step 1: estimation
C=zeros(T,N,N);
for t=1:T
Miss = isnan(X(t,:));
Obs = ~Miss;
c=zeros(N,N);
y = X(t,:)';
if ~isempty(X(Miss))
y(Miss) = M(Miss)+S(Miss,Obs)*inv(S(Obs,Obs))*(y(Obs)-M(Obs));
c(Miss,Miss) = S(Miss,Miss)-S(Miss,Obs)*inv(S(Obs,Obs))*S(Obs,Miss);
end
Y(t,:) = y';
C(t,:,:) = c + (y-M)*(y-M)';
end

% Step 2: update
M_new = mean(Y)';
S_new = squeeze(mean(C,1));

D4=[(M_new - M).^4;
diag((S_new - S).^2) ];
Distance= mean(D4.^(1/4));
Convergence = (Distance < Tolerance);

M=M_new;
S=S_new;
end

E_EM=M;
S_EM=S;```