python code(numpy torch)
import numpy as np

# a. answer
# Generate 100,000 consumer random disturbance items according to the first-type extreme value distribution
ee = np.random.gumbel(loc=0.0, scale=1.0, size=(100000, 4))
fake_data = np.zeros((100000, 4))
car1 = [1.6, 5.0]
car2 = [1.4, 3.0]
car3 = [1.5, 2.0]

# Generate the fake data
fake_data[:, 0] = ee[:, 0]
fake_data[:, 1] = -1.0 + 0.5 * car1[0] - 0.1 * car1[1] + ee[:, 1]
fake_data[:, 2] = -1.0 + 0.5 * car2[0] - 0.1 * car2[1] + ee[:, 2]
fake_data[:, 3] = -1.0 + 0.5 * car3[0] - 0.1 * car3[1] + ee[:, 3]

# b. answer
# Calculate market share
share1 = np.sum(np.argmax(fake_data, axis=1) == 1) / 100000
share2 = np.sum(np.argmax(fake_data, axis=1) == 2) / 100000
share3 = np.sum(np.argmax(fake_data, axis=1) == 3) / 100000

# c. answer
# Estimate the model useing the fake data
from sklearn.linear_model import LinearRegression
model = LinearRegression()
y = (np.mean(fake_data, axis=0) - np.mean(ee, axis=0))[1:]
model.fit(X=np.array([[1.6, 5.0], [1.4, 3.0], [1.5, 2.0]]), y=y)
a = model.intercept_
b = model.coef_
beta0 = a
beta1 = b[0]
beta2 = b[1]
print(beta0)
print(beta1)
print(beta2)

# d. answer
# Calculate the information matrix (Hessian matrix)
import torch
def get_jacobian(fun, x, noutputs):
    x_size = list(x.size())
    x = x.unsqueeze(0).repeat(noutputs, *([1]*(len(x_size))) ).detach().requires_grad_(True)
    y = fun(x)
    y.backward(torch.eye(noutputs))
    return x.grad.view(noutputs,*x_size)

def Hessian_matrix(fun, x):
    # y is a scalar Tensor
    def get_grad(xx):
        y = fun(xx)
        grad, = torch.autograd.grad(y, xx, create_graph=True, grad_outputs=torch.ones_like(y))
        return grad

    x_size = x.numel()
    return get_jacobian(get_grad, x, x_size)

weight = torch.rand(3, 2)

def fun(x):
    return (x @ weight).pow(2).sum(-1)

x = torch.tensor([beta0, beta1, beta2], requires_grad=True).float()
# the information matrix (Hessian matrix)
info = Hessian_matrix(fun, x)
print(info)
# compute parameter standard error
error = np.diag(np.linalg.inv(info.detach().numpy())) ** 0.5
print(error)

# e. answer
# Generate 100,000 consumer random disturbance items according to the first-type extreme value distribution
ee = np.random.gumbel(loc=0.0, scale=1.0, size=(100000, 4))
fake_data = np.zeros((100000, 4))
car1 = [1.6, 5.0]
car2 = [1.4, 3.0]
car3 = [1.5, 2.0]

# Generate the fake data
fake_data[:, 0] = ee[:, 0]
fake_data[:, 1] = -1.0 + 0.5 * car1[0] + np.random.lognormal(-0.1, 0.1, 1) * car1[1] + ee[:, 1]
fake_data[:, 2] = -1.0 + 0.5 * car2[0] + np.random.lognormal(-0.1, 0.1, 1) * car2[1] + ee[:, 2]
fake_data[:, 3] = -1.0 + 0.5 * car3[0] + np.random.lognormal(-0.1, 0.1, 1) * car3[1] + ee[:, 3]

from sklearn.linear_model import LinearRegression
model = LinearRegression()
y = (np.mean(fake_data, axis=0) - np.mean(ee, axis=0))[1:]
model.fit(X=np.array([[1.6, 5.0], [1.4, 3.0], [1.5, 2.0]]), y=y)
a = model.intercept_
b = model.coef_
beta0 = a
beta1 = b[0]
beta2 = b[1]
print(beta0)
print(beta1)
print(beta2)