function D=npest_instadv(y,z,x,i1,i2)

%function for the nonparametric estimation, simple program, z and x are
%endogenous and y and w in the second equation
%%%h1(y)=phi(z)+x+u
%%%h2(z)=psi(y)+w+v

n=length(y);

%bandwidths
hi1=1.59*std(i1)/(n^(1/5));
hi2=1.59*std(i2)/(n^(1/5));
hy=1.59*std(y)/(n^(1/5));
hz=1.59*std(z)/(n^(1/5));


dummy = (repmat(i1,1,n) - repmat(i1',n,1)) / hi1; Ki = kernel(dummy);
dummy = (repmat(i2,1,n) - repmat(i2',n,1)) / hi2; Ki = Ki .* kernel(dummy);
Ki = Ki ./ repmat(sum(Ki,2),1,n);

dummy = (repmat(z,1,n) - repmat(z',n,1)) / hz; Kzi = kernel(dummy);
dummy = (repmat(y,1,n) - repmat(y',n,1)) / hy; Kzi = Kzi .* kernel(dummy);
Kzi = Kzi ./ repmat(sum(Kzi,2),1,n);

Kz = (repmat(z,1,n) - repmat(z',n,1)) / hz; Kz =  kernel(Kz);
Kz = Kz ./ repmat(sum(Kz,2),1,n);

Ky = (repmat(y,1,n) - repmat(y',n,1)) / hy; Ky =  kernel(Ky);
Ky = Ky ./ repmat(sum(Ky,2),1,n);

%pmatrix
P=eye(n)-((1/n)*ones(n,n));

%definition of alpha parameters
Q=100;
alpha = logspace(-4,0,Q);  %alpha=alpha';
N_alpha = length(alpha);

for i_alpha=1:N_alpha;
  BMi=alpha(i_alpha)*eye(n)+Kzi*Ki ;
  sonuci=(Kzi)*Ki*x;
  estimate_1 = inv(BMi)*sonuci;
  estimate_2=inv(BMi)*(sonuci+estimate_1*alpha(i_alpha));
  res=sonuci-Kzi*Ki*estimate_2;
  Ea2(i_alpha) = norm(res)^2;
end;
critere=Ea2./(alpha.^2);
[y1,I1]=min(critere);
best_alpha11 = alpha(I1)
% best_alpha11 = 0.0739;
% % figure(1); loglog(alpha,critere,'.-', [best_alpha11],[critere(I1)],'rx');hold on;  
Q=100;
c = logspace(0,5,Q);  %alpha=alpha';
N_c = length(c);

sonuc=[Ky*Ki*x; P*Kz*Ki*x];
for i_c=1:N_c;
  BMc=[best_alpha11*eye(n)+Ky*Ki -Ky*Ki; P*Kz*Ki -((best_alpha11)*c(i_c)*eye(n)+P*Kz*Ki)]; ;
  estimate_c1 = inv(BMc)*sonuc;
  estimate_c2=inv(BMc)*(sonuc+[estimate_c1(1:n);estimate_c1(n+1:2*n)*c(i_c)]);
  resc=sonuc-([Ky*Ki -Ky*Ki; Kz*Ki -Kz*Ki]*estimate_c2);
  
  Eac(i_c) = norm(resc)^2;
end;

c11=c.*c;

criterec=Eac./c11;
[yc,Ic]=min(Eac);
best_c = c(Ic)
% best_c=1;
% figure(11); loglog(c,Eac,'.-', [best_c],[Eac(Ic)],'rx');hold on; 

BMb2=[best_alpha11*eye(n)+Ky*Ki -Ky*Ki; P*Kz*Ki -((best_alpha11*best_c)*eye(n)+P*Kz*Ki)];

% sonuc=[Ky*x; P*Kz*x];
best_estimate2 = inv(BMb2)*sonuc;


hy2(1:n)=best_estimate2(1:n);
phiden2(1:n)=best_estimate2(n+1:2*n);
D=[hy2' phiden2'];
