s=randn('state');

T=cumsum(randn(500,1)); % Random walk temperature (if not ok, what assumption that prevents this?)
R=randn(500,6)*10; % Random noise on proxies (if too much, what is the assumption for noise level?) 
Tm=T;

% 'RegEM assumes MAR, not MCAR' 
% 'All three cases are MAR'
% Trying cases 2 and 3:


Tm(101:400,:)=NaN;save rstate_middle s; % missing values at the middle (2), or
% Tm(201:500,:)=NaN;save rstate_end s;  % missing values at the end (3). I didn't check the values itself
 
 P=(repmat(T,1,6)+R)*diag(randn(6,1))*10; % pseudoproxies with random scale

%Y= incomplete data matrix, which contains both instrumental
%data (surface temperature grid-box values
%arranged with rows representing the years and columns
%representing grid boxes) and proxy data 
Y=[Tm P]; 

[Xe, M, C, Xerr] = regem(Y); % run the original regem without options

close all
plot(Xe(:,1),'g')
hold
plot(T,'r')
plot(Xe(:,1)+2*Xerr(:,1)) % Xerr= an estimated standard error of the imputed values 
plot(Xe(:,1)-2*Xerr(:,1))

legend('Imputed','True','+2se','-2se')