%% CCE asymptotic bias test
% see e.g. Sristava and Singh (1989) Small-Disturbance Asymptotic Theory
% for Linear-Calibration Estimators. Technometrics, August 1989, Vol.
% 31, NO. 3
si=30000; % how many simulations
Rc=zeros(si,1); % out
Rcvm=Rc;
A=100; % total number of samples
n=18; % calibration data 1:n
b=2; % unknown scale
x=(1:A)'; % x(n+1:end) will be estimated %
stdn=5; % noise std
for i=1:si
N=randn(A,1)*stdn; % i.i.d Gaussian noise
y=b*x+N; % 'proxy'
xc=x(1:n)-mean(x(1:n)); % centered reference data
yc=y(1:n)-mean(y(1:n)); % centered observations
Xc=mean(x(1:n))+(xc'*xc)/(xc'*yc)*(y-mean(y(1:n))); %% classical calibration estimator
Xcvm=mean(x(1:n))+sqrt(xc'*xc)/sqrt(yc'*yc)*(y-mean(y(1:n))); % CVM (my interpretation)
Resc=Xc-x; % estimation error
Rescvm=Xcvm-x;
Rc(i)=Resc(end); % store the last value
Rcvm(i)=Rescvm(end);
end
close all
figure
plot(Rc)
hold on
plot(Rcvm,'r')
legend('Res CCE','Res CVM')
disp('CCE Aymptotic bias:')
U=-stdn^2*(mean(x(1:n))-x(end));
D=(b*b*(xc'*xc));
U/D % from the reference, note that 'small as possible errors' are assumed
disp('CCE Sample bias:')
mean(Rc)
disp('CVM Sample Bias:')
mean(Rcvm)
figure
subplot(211)
hist(Rc,si/100)
AA=get(gca,'XLim');
subplot(212)
hist(Rcvm,si/100)
set(gca,'XLim',AA)
xlabel('CCE (top), CVM (bottom)')