# What is the equivalent Julia code for statsmodels.api.OLS

**URL:** <https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070>\
**Category:** Statistics\
**Tags:** python\
**Created:** [March 25, 2024, 9:16am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070 "2024-03-25T09:16:20Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 25, 2024, 9:16am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/1 "2024-03-25T09:16:20Z")

</div>

I want to write the equivalent Julia code for the following Python code:

```python
import statsmodels.api as sm
reg = sm.OLS(Y,sm.add_constant(X))
res = reg.fit()
res.params[1]/res.bse[1]

```

The copilot asserts that the equivalent Julia code is:

```julia
using GLM
reg = glm(X, Y, Normal(), IdentityLink())
coef(reg)[2]

```

However, the result of these two is not equal. Example:

```python
# Python
import numpy as np
import statsmodels.api as sm
X = np.array([[1, 2, 3],[4, 5, 6], [7, 8, 9]])
y = np.array([1, 0, 2])
reg = sm.OLS(y,sm.add_constant(X))
res = reg.fit()
res.params[1]/res.bse[1]
# -3.5606350123330477e-16

```

```julia
# Julia
using GLM
X = [
1 2 3
4 5 6
7 8 9
]
y = [1, 0, 2]
reg = glm(X, y, Normal(), IdentityLink())
coef(reg)[2]
# 0.0

```

Note that I am using `statsmodels` v0.13.5 in this case.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 25, 2024, 9:58am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/2 "2024-03-25T09:58:14Z")

</div>

Your system is singular, the coefficients are not estimable.

```julia
>>> res.summary()
...
...

[2] The input rank is higher than the number of observations.
[3] The smallest eigenvalue is 6.54e-32. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.

```

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 25, 2024, 10:11am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/3 "2024-03-25T10:11:07Z")

</div>

And, also, you are adding a constant to `X` in the python regression. You should do that in julia too. In your example it can by done by `X = hcat([1,1,1], X)`. If you use the formula interface (with `DataFrames` or another `Table`) rather than the matrix interface, it’s simpler.

```julia
using DataFrames
df = DataFrame(X, :auto)
df.y = y
fit(LinearModel, @formula(y ~ x1 + x2 + x3 + 1), df)

```

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 25, 2024, 11:06am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/4 "2024-03-25T11:06:39Z")

</div>

Just to be complete. Python uses a pseudo inverse for degenerate systems. There is really no canonical way to solve singular problems. Julia typically does it by setting superfluous coefficients to zero, with standard error `Inf` or `NaN`. You can mimic python’s behaviour as follows. However, with degenerate systems the interpretation of coefficients and standard errors will anyway be wishy-washy.

```julia
using LinearAlgebra

Xc = hcat([1,1,1], X) # add a constant
cvmat = pinv(Xc' * Xc) # compute covariance matrix as pseudoinverse
coefs = cvmat * Xc' * y # perform the OLS
serr = sqrt.(diag(cvmat) .* 3/2) # find standard errors, small sample adjust by N/(N-dof)
coefs[2]/serr[2] # compute t-score

```

---

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 25, 2024, 11:35am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/5 "2024-03-25T11:35:10Z")

</div>

> [@sgaure](#):
>
> Just to be complete. Python uses a pseudo inverse for degenerate systems. There is really no canonical way to solve singular problems. You can mimic python’s behaviour as follows:
> 
> ```julia
> using LinearAlgebra
> 
> Xc = hcat([1,1,1], X) # add a constant
> cvmat = pinv(Xc' * Xc) # compute covariance matrix as pseudoinverse
> coefs = cvmat * Xc' * y # perform the OLS
> serr = sqrt.(diag(cvmat) .* 3/2) # find standard errors, small sample adjust by N/(N-dof)
> coefs[2]/serr[2] # compute t-score
> 
> ```

Thank you so much for your replies. I don’t understand why you’ve used the constant `3/2`. Anyway, the result is still different:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/9/59423303866328f4c63ceef8b5bd3c6db4027678.png)

Am I missing something here?

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 25, 2024, 11:44am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/6 "2024-03-25T11:44:30Z")

</div>

> [@Shayan](#):
>
> Am I missing something here?

Yes. You get something very close to zero. Something like 0.00000000000000002. Both python’s and julia’s numbers are essentially equal to zero. If you look at the other coefficients and standard errors, they are essentially equal.

Btw, the multiplication by 3/2 is a small sample correction of the standard errors. I.e. `N/(N-1)` where `N` is the number of observations. It does normally not have any large impact since the number of observations `N` is typically in the hundreds (or thousands, or millions). But in your example `N==3`.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 25, 2024, 1:21pm UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/7 "2024-03-25T13:21:31Z")

</div>

Here’s an interesting thing with the python script:

```julia
import numpy as np
import statsmodels.api as sm
X = np.array([[1, 2, 3],[4, 5, 6], [7, 8, 9]])
y = np.array([1, 0, 2])
reg = sm.OLS(y,sm.add_constant(X))
res = reg.fit()
res.params[1]/res.bse[1] # this is some small number
res = reg.fit(method='qr') # this one fails 
res = reg.fit() # do it again with default method 'pinv'
res.params[1]/res.bse[1]
# returns nan (because res.bse[1] is nan as it should be)

```

---

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 26, 2024, 10:32am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/8 "2024-03-26T10:32:06Z")

</div>

> [@sgaure](#):
>
> You get something very close to zero. Something like 0.00000000000000002. Both python’s and julia’s numbers are essentially equal to zero.

But the signs are different. The output in the Julia is a positive value while the Python’s output is a negative value. There should be something different for sure. So, I don’t think the proposed Julia code is equivalent to the Python implementation.

> [@sgaure](#):
>
> If you look at the other coefficients and standard errors, they are essentially equal.

Let me check this out. If the coefficients are similar, then there should be the equivalent Julia code for the `res.params[1]/res.bse[1]`. Otherwise, the problem lies in fitting the OLS in Julia implementation in the same way as `sm.OLS` does. Actually, this `res.params[1]/res.bse[1]` is the important part of the code to me and the final result that should be used in my project.

---

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 26, 2024, 10:35am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/9 "2024-03-26T10:35:22Z")

</div>

> [@sgaure](#):
>
> `returns nan (because res.bse[1] is nan as it should be)`

I didn’t get the point. What is the conclusion?

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 26, 2024, 11:58am UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/10 "2024-03-26T11:58:43Z")

</div>

The main point is that OLS on a singular system is not well defined. As I have told you numerous times. The coefficients are not estimable. Entirely arbitrary choices must be made to get some coefficients. I have described julia’s choices, and python’s choices, and described how you can mimic python’s choices in julia. The coefficients you obtain are relative to these choices, and typically are not meaningful without analyzing the nature of the singularity (or multicollinearity which it is often called in this context). Though the model as such can be used for prediction.

And, of course, as always with floating point operations, when comparing results it must be done up to a precision, there are no meaningful exact comparisons of floating point numbers, except in some very simple cases. The problem at hand, computing a pseudo inverse, involves inverting singular values which are not too small, i.e. one must choose what “too small” means, this will have an impact on the pseudoinverse, and consequently on the computed coefficients. Julia and python may have slightly different parameters for computing the pseudoinverse, and even slightly different ways to compute the singular value decomposition. I don’t know.

Also, in the python-mimicry I multiplied the covariance matrix with `Xc' * y` to find the coefficients. Because I needed the covariance matrix anyway to find the standard errors. It is more common to solve the system `Xc' * Xc * x = Xc' * y` for `x`, which is mathematically the same, but typically uses a different method, e.g. a pivoted Cholesky decomposition or qr or whatnot. This will yield slightly different estimates, i.e. in the 17th digit or so. What exactly python does, I don’t know. Fortunately, it doesn’t matter, because of the following.

If your program relies on exact floating point values which varies with more or less arbitrary choices made by some package, you are doing it the wrong way.

Now, what you are computing is known as “t-score” or “t-value”, how many standard errors your estimate is away from zero. This is a sensible thing to compute, and is used to judge whether the results are statistically significant. When the t-score is small (how small depends on your application), the coefficient is indistinguishable from zero with the data you have. The difference between `0.001` and `-0.001`, with a standard error around `1` is simply completely irrelevant.

---

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 26, 2024, 6:10pm UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/11 "2024-03-26T18:10:17Z")

</div>

Ok, this is going to be a long reply:  
The aim is not to interpret the t-score! The value itself matters. I want to translate a Python code that detects bubble prices in historical data of an asset (like Microsoft, AKA MSFT). **The code I was asking for is a part of this procedure**. I will provide the full Python code and the translated Julia code in the following section to show why the value of `res.params[1]/res.bse[1]` matters not the interpretation of it. I run both codes on the same data to compare the results. I compared the `X` and `y` in both scripts to make sure the same data is being fed to the OLS model. I did it by printing the values of `X` and `y` and raising an arbitrary error after the print statements (to stop the run and compare the results by observing them).

> **Pyhton and Julia implementation in a glance**
>
> ![image](https://global.discourse-cdn.com/julialang/original/3X/3/e/3e551d715130100b3a690d844e04c21fdb5b48f1.jpeg)

> **Pyhton implementation**
>
> ```python
> import numpy as np
> import pandas as pd
> import statsmodels.api as sm
> import matplotlib.pyplot as plt
> prices = pd.read_csv("../data/GDP.csv", header=0).GDP
> r0 = int(len(prices)*0.1)
> #specify lags for the augmented Dickey-Fuller (ADF) test
> adf_lags = 3
> #critical value of the right-tailed ADF-test (95%) from Phillips et al. (2015)
> crit = 1.49
> #transforming data
> log_prices = np.array(np.log(prices))
> delta_log_prices = log_prices[1:] - log_prices[:-1]
> n = len(delta_log_prices)
> BSADF = np.array([])
> #calculating ADF stats
> for r2 in range(r0,n):
> ADFS = np.array([])
> for r1 in range(0,r2-r0+1):
> X0 = log_prices[r1:r2+1]
> X = pd.DataFrame()
> X[0] = X0
> for j in range(1,adf_lags+1):
> X[j] = np.append(np.zeros(j),delta_log_prices[r1:r2+1-j])
> X = np.array(X)
> Y = delta_log_prices[r1:r2+1]
> reg = sm.OLS(Y,sm.add_constant(X))
> res = reg.fit()
> ADFS = np.append(ADFS, res.params[1]/res.bse[1])
>         
> BSADF = np.append(BSADF, max(ADFS))
> #visualising the results
> plt.rc('xtick',labelsize = 8)
> plt.plot(prices.index[r0+1:],BSADF)
> plt.plot(prices.index[r0+1:],np.ones(len(BSADF))*crit)
> #printing dates when bubbles were detected
> print(prices.index[r0+1:][BSADF > crit])
> 
> ```

> **Julia implementation**
>
> ```julia
> using LinearAlgebra
> using Plots
> 
> function bubbledetector(lnprices::AbstractVector, adf_lags::Integer, crit::AbstractFloat)
> Δlnprices = diff(lnprices)
> n = length(Δlnprices)
> r0 = round(Int, length(lnprices)*0.1)
> BSADF = similar(lnprices, n-1-r0+1)
> for (idxr2, r2) ∈ enumerate(r0:n-1)
> ADFS = similar(lnprices, r2-r0+1)
> for (idxr1, r1) ∈ enumerate(0:r2-r0)
> X₀ = lnprices[r1+1:r2]
> X = similar(lnprices, length(X₀), adf_lags+1)
> X[:, 1] = X₀
> for j in 1:adf_lags
> X[:, j+1] = vcat(zeros(j), Δlnprices[r1+1:r2-j])
> end
> Y = Δlnprices[r1+1:r2]
> Xc = hcat(ones(size(X, 1)), X) # add a constant
> cvmat = pinv(Xc' * Xc) # compute covariance matrix as pseudoinverse
> coefs = cvmat * Xc' * Y # perform the OLS
> serr = sqrt.(diag(cvmat) .* length(Y)/(length(Y)-1)) # find standard errors, small sample adjust by N/(N-dof)
> ADFS[idxr1] = coefs[2]/serr[2]
> end
> BSADF[idxr2] = maximum(ADFS)
> end
> return BSADF, findall(BSADF .> crit)
> end
> 
> function plotbubles(BSADF::AbstractVector, idxs::AbstractVector, crit::AbstractFloat)
> plot(BSADF, label="BSADF")
> plot!(ones(length(BSADF))*crit, color=:red, linestyle=:dash, label="Critical Value")
> plot!(idxs, BSADF[idxs], seriestype=:scatter, color=:red, label="Bubbles", alpha=0.5)
> end
> 
> prices = identity.(readdlm("./data/GDP.csv", ',', header=true)[1][:, 2])
> BSADF, idxs = bubbledetector(log.(prices), 3, 0.8)
> plotbubles(BSADF, idxs, 0.8)
> 
> ```

> **Python's and Julia's output in a glance**
>
> ![image](https://global.discourse-cdn.com/julialang/original/3X/3/9/39412d7be11508bd0587947c2b53658eb97fdaf2.png)  
> The thing is, the values/trend isn’t even similar to the Python’s output! Not even close! Meaning the values are extremely different! As I said earlier, I checked the inputs (`X` and `y`) in both scripts to make sure the same data is being fed to the OLS models. In conclusion, the difference is because of the OLS implementation.

I brought up all these details to explain why I need the same output as `res.params[1]/res.bse[1]`.

---

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 26, 2024, 6:20pm UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/12 "2024-03-26T18:20:24Z")

</div>

> [@sgaure](#):
>
> The coefficients are not estimable. Entirely arbitrary choices must be made to get some coefficients.

So, there is no hope to get the same result necessarily due to different implementations. Thank you.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 26, 2024, 7:00pm UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/13 "2024-03-26T19:00:31Z")

</div>

There is no hope of getting meaningful results with a singular system such as in your example. It seems your real system is not degenerate, and then you can of course get the same results, up to the hardware precision. You can use the `fit(LinearModel, ...)` as I suggested, or `glm(...)` as you tried, but with the intercept added, i.e. `glm(Xc, y, Normal(), IdentityLink())`.

I thought your main issue was that you wanted results down to the same bit when the design matrix is singular, as in your example. That’s a tall order.

Now, anyway, I see I made an error in my julia program. Python does not seem to do a small sample correction of the standard errors, so that one can be removed. By a coincidence in the example that factor happened to be equal to the one I forgot, to multiply by the residual variance, so the standard errors can be computed with something like this:

```julia
resid = norm(Xc * coefs .- y)
serr = sqrt.(diag(cvmat)) .* resid

```

And, actually, if you only use a single element, you can compute it as `serr2 = sqrt(cvmat[2,2]) * resid`.

---

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 26, 2024, 8:14pm UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/14 "2024-03-26T20:14:41Z")

</div>

> [@sgaure](#):
>
> There is no hope of getting meaningful results with a singular system such as in your example. It seems your real system is not degenerate, and then you can of course get the same results, up to the hardware precision.

Yes, I got your point.

> [@sgaure](#):
>
> I thought your main issue was that you wanted results down to the same bit when the design matrix is singular, as in your example. That’s a tall order.

You’re right. I should have brought up the details in the first place. I am sorry for any confusion or inconvenience that may caused.

> [@sgaure](#):
>
> You can use the `fit(LinearModel, ...)` as I suggested, or `glm(...)` as you tried, but with the intercept added, i.e. `glm(Xc, y, Normal(), IdentityLink())`.

Considering `reg = glm(Xc, Y, Normal(), IdentityLink())`, I don’t know what should be the next line to compute the equivalent `res.params[1]/res.bse[1]` in this case.

> [@sgaure](#):
>
> And, actually, if you only use a single element, you can compute it as `serr2 = sqrt(cvmat[2,2]) * resid`.

This is great. Although there are some negligible differences as shown in the picture below, the output’s trend is extremely similar to Python’s output. The only thing is the range of output. It seems like in Python, the values are multiplied by a constant value \> 1 (it seems that the final `BSADF` values in Julia should be multiplied by 10):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/b/3b68fc20db97fe820b289a59ea1e6cf87cc35a27.png)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/1/614445320b674cd177c2edbf55e743c9529cb888.png)

> **The following Julia code is used for the plot above**
>
> ```julia
> using LinearAlgebra
> using Plots
> 
> function bubbledetector(lnprices::AbstractVector, adf_lags::Integer, crit::AbstractFloat)
> Δlnprices = diff(lnprices)
> n = length(Δlnprices)
> r0 = round(Int, length(lnprices)*0.1)
> BSADF = similar(lnprices, n-1-r0+1)
> for (idxr2, r2) ∈ enumerate(r0:n-1)
> ADFS = similar(lnprices, r2-r0+1)
> for (idxr1, r1) ∈ enumerate(0:r2-r0)
> X₀ = lnprices[r1+1:r2]
> X = similar(lnprices, length(X₀), adf_lags+1)
> X[:, 1] = X₀
> for j in 1:adf_lags
> X[:, j+1] = vcat(zeros(j), Δlnprices[r1+1:r2-j])
> end
> Y = Δlnprices[r1+1:r2]
> Xc = hcat(ones(size(X, 1)), X) # add a constant
> cvmat = pinv(Xc' * Xc) # compute covariance matrix as pseudoinverse
> coefs = cvmat * Xc' * Y # perform the OLS
> resid = norm(Xc * coefs .- Y)
> serr = sqrt(cvmat[2,2]) * resid
> ADFS[idxr1] = coefs[2]/serr
> end
> BSADF[idxr2] = maximum(ADFS)
> end
> return BSADF, findall(BSADF .> crit)
> end
> 
> function plotbubles(BSADF::AbstractVector, idxs::AbstractVector, crit::AbstractFloat)
> plot(BSADF, label="BSADF")
> plot!(ones(length(BSADF))*crit, color=:red, linestyle=:dash, label="Critical Value")
> plot!(legend=:outerright)
> end
> 
> prices = identity.(readdlm("./data/GDP.csv", ',', header=true)[1][:, 2])
> BSADF, idxs = bubbledetector(log.(prices), 3, 1.49)
> plotbubles(BSADF, idxs, 1.49)
> 
> ```

In summary, I used the following lines as the equivalent Julia code for `sm.OLS`:

```julia
Xc = hcat(ones(size(X, 1)), X) # add a constant
cvmat = pinv(Xc' * Xc) # compute covariance matrix as pseudoinverse
coefs = cvmat * Xc' * Y # perform the OLS
resid = norm(Xc * coefs .- Y)
serr = sqrt(cvmat[2,2]) * resid
coefs[2]/serr

```

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [March 26, 2024, 10:01pm UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/15 "2024-03-26T22:01:04Z")

</div>

> [@Shayan](#):
>
> Considering `reg = glm(Xc, Y, Normal(), IdentityLink())`, I don’t know what should be the next line to compute the equivalent `res.params[1]/res.bse[1]` in this case.

It’s here: [Home · GLM](https://juliastats.org/GLM.jl/stable/#Methods-applied-to-fitted-models)

```julia
coef(reg)[2] / stderror(reg)[2]

```

> The only thing is the range of output. It seems like in Python, the values are multiplied by a constant value \> 1 (it seems that the final `BSADF` values in Julia should be multiplied by 10):

This I do not understand. However, I took the OLS equations from the top of my head, and they work fine in a toy example:

```julia
using LinearAlgebra, GLM
X = [
    1 3 4
    2 8 1
    4 5 2
    6 7 2
    9 4 3
]

y = [5,4,6,2,3]
Xc = hcat(ones(size(X,1)), X)

# GLM:

regglm = glm(Xc, y, Normal(), IdentityLink())

# easy fit

reglm = fit(LinearModel, Xc, y)

#manual
cvmat = inv(Xc' * Xc)
cf = cvmat * Xc' * y
resid = norm(Xc * cf .- y)
se = sqrt.(diag(cvmat)) * resid

println("GLM:\n", regglm)
println("\n\nLM:\n", reglm)
println("\n\nMAN:")
[cf se]

```

and the same in python:

```julia
import numpy as np
import statsmodels.api as sm

X = np.array([[1, 3, 4],[2,8,1], [4,5,2], [6,7,2], [9,4,3]])
y = np.array([5,4,6,2,3])
reg = sm.OLS(y,sm.add_constant(X))
res = reg.fit()
res.params
res.bse

```

They are all equal.

---

<div class="post-metadata">

**Author:** ![Shayan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shayan/32/33025_2.png) [@Shayan](https://discourse.julialang.org/u/Shayan)\
**Post date:** [March 26, 2024, 10:17pm UTC](https://discourse.julialang.org/t/what-is-the-equivalent-julia-code-for-statsmodels-api-ols/112070/16 "2024-03-26T22:17:58Z")

</div>

> [@sgaure](#):
>
> `coef(reg)[2] / stderror(reg)[2]`

Thank you so much. This is exactly what I needed and the output exactly mimics the `sm.OLS`’s output! Thank you so much for your help.