function A = aoyagi11func(B, maxa, mina, ia, na, M, delta, edm, mu)

nb = 4; %number of buyers

% define grid for z
maxz = 1; 						% maximum value of z grid
minz = 0; 						% minimum value of z grid
iz = 0.005; 					% size of z grid
nz = round((maxz - minz)/iz+1); % number of z grid points

% define z as a column vector
for i=1:nz
	z(i,1) = (i-1)*iz + minz;
end

% set parameter values
c0 = 0;
ci = ones(nb,1);
for i = 1 : nb;
	ci(i,1) = B(1,i);
end;

%initialize payoff vectors
pi = ones(nz,na); 		% values when (z,a)
pistar = ones(nb,na); 	% optimal value when state is a
pistara = ones(nz,1);	% expected optimal value when accepted
pistarr = ones(nz,1);	% expected optimal value when rejected

%t=4 start from the final period
t=4;
maxat = (t-1)*(M-mu);	% maximum of alpha in t
minat = (t-1)*(-mu);	% minimum of alpha in t
for i = 1:na
	a = mina + (i-1)*ia;		% compute alpha for the ith grid point

	for j = 1:nz;
		
		theta = -1/delta * log(1-(1-edm)*(1-z(j,1))); % theta for each z
		edth = exp(-delta*theta);
		
		if a > maxat | a < minat	% if the ith grid point exceeds the bounds
			pi(:,i) = -100*nb;		% substitute a small value
		else
			pi(j,i) = z(j,1)*theta + ci(t)*a*z(j,1) + c0*z(j,1);
		end
		
	end;
	
end

pistar(t,:) = max(pi);

%t=1,2,3
t=3;
while t>0
	maxat = (t-1)*(M-mu);
	minat = (t-1)*(-mu);
	for i = 1:na
		a = mina + (i-1)*ia;
		
		k = ones(nz,1); % kappa(z) as a nz vector
		l = ones(nz,1); % lambda(z) as a nz vector
		
		for j = 1:nz;
		
			theta = -1/delta * log(1-(1-edm)*(1-z(j,1))); % theta for each z
			edth = exp(-delta*theta);
			
			if z(j,1) == 0
				k(j,1) = M-mu;
			else
				k(j,1) = ((theta + delta^(-1))*edth-(M+delta^(-1))*edm) / (edth-edm)...
					- (delta^(-1)-(M+delta^(-1))*edm) / (1-edm);
			end;
			
			if z(j,1) == 1
				l(j,1) = -mu;
			else
				l(j,1) = (delta^(-1)-(theta+delta^(-1))*edth) / (1-edth)...
					- (delta^(-1)-(M+delta^(-1))*edm) / (1-edm);
			end;
			
			aa(j) = a+k(j,1);	% future state for each z when accepted
			ar(j) = a+l(j,1);	% future stete for each z when refected
			
			if aa(j) > t*(M-mu);
				aa(j) = t*(M-mu);
			end;
			if ar(j) < t*(-mu);
				ar(j) = t*(-mu);
			end;
			
			gaa = fix((aa(j) - mina)/ia) + 1; % grid point for aa(j)
			gar = fix((ar(j) - mina)/ia) + 1; % grid point for ar(j)
			
			pistara(j,1) = pistar(t+1,gaa); % future payoff pistar_a when accepted
			pistarr(j,1) = pistar(t+1,gar); % future payoff pistar_r when rejected
			
			if a > (maxat+ia) | a < (minat-ia)
				pi(j,i) = -100*nb;
			else
				pi(j,i) = z(j,1)*theta + (ci(t)*a+c0)*z(j,1) + ...
						z(j,1)*pistara(j,1) + (1-z(j,1))*pistarr(j,1);
			end
			
		end
	end
	
	pistar(t,:) = max(pi); % choose the maximizing z for each a in period t
	t=t-1;
	
end
A(1,:) = pistar(1,:);

