1. Write a function called sys_internal(x) that expects a 2x1 vector and returns $g(x).$ Test your code by sending the vector $\l(\matrix{1\cr 1}\r)$ to sys_internal() and printing the output. You should get g():

         -1.0000
         0.30117
   

2. Write a new function named sys() that follows the format required by SolveNLE and is a "wrapper" around your simpler sys_internal() function:

   sys(ag,x)

   It returns either 0 if $x_0<0$, otherwise it calls sys_internal() and returns 1. Verify your sys computes the correct value for <1;1> and returns 1.
3. Use SolveNLE to solve for a root of the system. Check the convergence code. Try different starting values to see if it finds more than one root. Start at an undefined value (say $x_0=-1$) and see what happens.

1. Solve this problem numerically using SolveNLE() and Broyden's method following examples already done. Start with use of a numerical Jacobian.
2. Write your main code to start at values near the two solutions and see if they converge to the closest one.
3. Code the analytic Jacobian and run Newton-Raphson with it.
4. If you increase $v$ from 3 the two roots move closer together until eventually they meet when the line is tangent to the circle. If you go past that point there is no solution. Experiment with changing $v$ to find the tangency and see what happens when no solution exists. Here is a loop that would make it easier to do this ... insert in your main code and adjust $v$. You might also confirm analytically the value of $v$ that results in a unique root and enter it to confirm.

   while (TRUE) {
       scan("Enter v ","%g",&v);
       // reset your starting/ vector here
       SolveNLE( ...  );
       }
       

We can get back to the simple logit by setting J=2 and $\beta_0 = \overrightarrow{0}.$ That is option $0$ has beta vector of 0s, which means $e^{v_{i 0}} = e^{x_i\beta_0} = e^0 = 1.$ And now that function we've been working with all along becomes a special case of the multinomial logit log-likelihood function: $$L^J = \sum_{i=0}^{N-1} \ln \left( { e^{v_{i\,y_i}} \over \sum_{j=0}^{J-1} e^{v_{i\,j}} }\right),\qquad\hbox{and } v_{ij} \equiv x_i\beta_j.$$ Note that $L^J$ is shorthand. A more accurate definition would be $L^J(\beta_1,\dots,\beta_{J-1}; y, X).$ That is, $y$ and $X$ are taken as given. The vectors of coefficients are variable in order to maximize $L^J.$

In the simple logit we could maximize this function using Ox's MaxBFGS() by sending a single beta vector as the parameters to maximize over. In the multinomial model we want to maximize over 2 or more vectors, but MaxBFGS() like almost all optimizers takes only a column vector as the argument to maximize over. Let's define $\theta$ as that vector and think of it containing the "concatenation" of the $\beta_j$'s: $$\theta \equiv \pmatrix{\beta_1\cr\beta_2\cr\vdots\cr\beta_{J-1}}$$ Note that we don't include $\beta_0$ because it is fixed at 0 and MaxBFGS() should not change its values. A function written to compute $L^J$ would "unpack" $\theta$ into a matrix: $$\theta_{(J-1)K\times 1} \qquad \Rightarrow \qquad B_{K\times J} = \pmatrix{ \overrightarrow{0} & \beta_1 & \beta_2 & \cdots & \beta_{J-1}}.$$ Now notice that $XB$ will contain a $N\times J$ matrix consisting of all the values $x_i\beta_j$ for all $i$ and all $j$. So computing $v_{ij}$ and ultimately $L^J$ can be done with vector functions (or with loops over $j$ or both $i$ and $j$). The code and data file in the mlogit folder is designed to build a function to compute $L^J$ and finding the maximizing coefficients using one of the Max????() algorithms. The data come from the example used in Stata to illustrate its mlogit command. The data and estimates are summarized with these screen captures:

Add code to the files in the folder to:

1. Complete coding of LJ
2. Create a starting theta vector
3. Maximize LJ (use algorithm of your choice)
4. Print out results to compare with Stata output
5. EXTRA:
1. copy over stata_logit_data.dta from the logit folder
2. use of loaddemo to load the simple logit data data
3. Reset K and J (*NOTE: do not decrement y since it is already 0/1)
4. Verify LJ() to provides same results as simple logit

1. Simulate $Y_N$ in \eqref{YN} which is the scaled deviation from the population mean. Remember you can use built-in functions and operators to compute the mean of each column of a matrix. Here is the Ox graphing function that will do that if Y contains a $1\times R$ row matrix:

       #include 
          // Draw simulated sample
          // Compute values of YN in Y
       	DrawDensity(0, Y, "$Y_N$", 1, 1,1);
   	SaveDrawWindow("YN.pdf");
   

   The last 3 1's are optional but should make the graph look the way you want ("Graphic Reference" describes the Ox drawing routines).
2. Not everything we do with sample data follows a normal distribution in large samples even with IID sampling. One example is what are called extreme statistics (whereas the mean is a central statistic). For example, for a sample of $X_i$ of size $N$ define as the largest value in the sample: $X^{max}_N \equiv \max_i\ X_i$ Using the same simulated sample of $\chi^2_3$, graph $X^{max}_N$ and note whether it "looks normal" for large $N$

1. Use the loaddemo() code to read in the $X$ matrix from the Stata sample. Instead of using the 3 outcome $y$ vector from the data, simulated new $y$ vectors multiple times using Monte Carlo. To do this, you need to have the "true" or population parameter vector, often denoted with a "0," so here $\theta_0.$ Set $\theta_0$ equal to the MLE estimates from the actual data: $\theta_0 = \hat\theta^{mle}.$ Write a function to simulate a sample of $y$ using the same $X$ matrix and $\theta_0.$ THat is compute the simulated $u_{ij}$ and compute $j^\star_i.$ That means you have to store the index of the element in a row or column vector that holds the max (in this case the choice with the highest utility for observation $j$). Ox has a maxcindex() function to do this (or use loops). Note it does not provide a maxrindex() so you may have to transpose your matrix to use it.
2. For a simulated $y$ vector re-estimate the coefficients using MaxBFGS() as in the previous assignment. Then put this code into a loop to estimate from $R$ simulated samples.
3. You can put all the simulated results in a matrix that is $2K \times R.$ Each column is a simulated MLE vector. Compute the distribution of the estimates and compare them to $\theta_0$. (The moments() function is a good option for this.)