Anomally with nleqslv function in R

Question

I am using nleqslv function in R to solve 10 equations. However, sometimes it gets stuck when it can't resolve some of the inputs. Here is an example of the code, please note that the for loop is solving the equations and for this example, I have hard-coded the input for 1 row only.

The for loop is where the equations are being solved. I use a divisor to optimize the function's performance, but in this example, the divisor is ineffective. When I remove the divisor, the equations can be solved, but the results are incorrect. The purpose of the divisor is to ensure the equations resolve to the expected results.

library(nleqslv)
library(dplyr)

# df_pe1 <-  read.csv("C:\\Users\\IEC\\OneDrive - International Economics Consulting Ltd\\2023\\Partial Equilibrium\\PE test results to WITS SMART\\JG_inputs.csv")


df_pe1 <- data.frame(A= c("A"),
                     B = "B",
                     C = "888",
                     D = "XXX",
                     E = 73792,
                     F = 10.1,
                     G= 1,
                     H = 0.998,
                     I = -7.6004338)

names(df_pe1) <- c("Reporter","Partner","HS6","Commodity","imp_world","imports","initial_tariff","initial_ROW_tariff","imp_elasticity")



df1 <- df_pe1 #%>% filter(HS6 %in% c(80231))



df1 <- df1 %>%
  mutate(diviser = case_when( imports >= imp_world ~ 1*10^(nchar(round(imp_world,0) )+1),
                              imports < imp_world ~ 1*10^(nchar(round(imports,0) )+1)
                              ))


df1 <- df1 %>% 
  mutate(imp_elasticity = case_when(imp_elasticity <= -3 ~ -1.5,
                                    .default = imp_elasticity))




df1[,5] <-df1[,5]/df1$diviser
df1[,6] <-df1[,6]/df1$diviser


df1$ROW_PMxM=df1$imp_world
df1$ROW_TWA=df1$initial_ROW_tariff *100

df1$ROW_Tnew <-df1$initial_ROW_tariff *100

df1$PMxM=df1$imports
df1$TWA<-df1$initial_tariff *100
df1$Tnew <- 0

df1$MELAST =  df1$imp_elasticity #demand elas

df1$XELAST <- 20 #supply elas
df1$PM0<-1

df1$MSUB<- 3#Substitution  els 
#df1$substitution_elasticity

df1$HS<-df1$HS6

df1$M0<-df1$PMxM/df1$PM0
df1$ROW_M0<-df1$ROW_PMxM/df1$PM0

df1$PMA=1

df1$PX0<-df1$PM0/(1+df1$TWA/100)
df1$ROW_PX0<-df1$PM0/(1+df1$ROW_TWA/100)

df1$X0<-df1$M0
df1$ROW_X0<-df1$ROW_M0

df1$MA0<-df1$PMxM+df1$ROW_PMxM
df1$MSUBPAR<--1/df1$MSUB+1

df1$MSHIFT<-df1$MA0/df1$PMA^df1$MELAST


df1$XSHIFT<-df1$X0/df1$PX0^df1$XELAST
df1$ROW_XSHIFT<-df1$ROW_X0/df1$ROW_PX0^df1$XELAST

df1$Mshare<-(1/df1$M0^(df1$MSUBPAR-1))/(1/df1$M0^(df1$MSUBPAR-1)+1/df1$ROW_M0^(df1$MSUBPAR-1))
df1$ROW_Mshare<-(1/df1$ROW_M0^(df1$MSUBPAR-1))/(1/df1$M0^(df1$MSUBPAR-1)+1/df1$ROW_M0^(df1$MSUBPAR-1))

df1$Ashift<-df1$MA0/(df1$Mshare*df1$M0^df1$MSUBPAR+df1$ROW_Mshare*df1$ROW_M0^df1$MSUBPAR)^(1/df1$MSUBPAR)

is.nan.data.frame <- function(x)
  do.call(cbind, lapply(x, is.nan))

df1[is.nan(df1)] <- 0


dftemp<-df1



df<-dftemp

dftemp1<-df



dftemp2<-dftemp1 #1ST BATCH


results<-data.frame(matrix(nrow=nrow(dftemp2), ncol=10))
names(results)<- c("HS","X","PX","ROW_X","ROW_PX","M","ROW_M","PM","ROW_PM","PMA1")
names(dftemp2)
fn3 <- function(x) {
  eq1<-x[1]-df["XSHIFT","Partner"]*x[2]^df["XELAST","Partner"] #X_SUP(P) CIS
  eq2<-x[3]-df["XSHIFT","ROW"]*x[4]^df["XELAST","ROW"] #X_SUP(P) ROW
  eq3<-x[1]-x[5]                                       #EQUIL(P) CIS
  eq4<-x[3]-x[6]                                       #EQUIL(P) ROW
  eq5<-x[7]-x[2]*(1+df["T","Partner"]/100)                 #TARIFF(P) CIS
  eq6<-x[8]-x[4]*(1+df["T","ROW"]/100)                 #TARIFF(P) ROW
  
  eq7<-x[10]-df["ASHIFT",2]*(df["MSHARE","Partner"]*x[5]^df["MSUBPAR",2]+
                               df["MSHARE","ROW"]*x[6]^df["MSUBPAR",2])^(1/df["MSUBPAR",2]) #ARM (Sum ROW and CIS)
  
  eq8<-x[7]-x[9]*df["ASHIFT",2]*(df["MSHARE","Partner"]*(x[5]^df["MSUBPAR",2])+
                                   df["MSHARE","ROW"]*(x[6]^df["MSUBPAR",2]))^((1/df["MSUBPAR",2])-1)*
    df["MSHARE","Partner"]*(x[5]^(df["MSUBPAR",2]-1))                                                      #M_DEM(P) (Sum ROW and CIS) import demand by CIS
  
  eq9<-x[8]-x[9]*df["ASHIFT",2]*(df["MSHARE","Partner"]*(x[5]^df["MSUBPAR",2])+
                                   df["MSHARE","ROW"]*(x[6]^df["MSUBPAR",2]))^((1/df["MSUBPAR",2])-1)*
    df["MSHARE","ROW"]*(x[6]^(df["MSUBPAR",2]-1))                                                      #M_DEM(P) (Sum ROW and CIS) import demand by ROW
  
  eq10<-x[10]- df["MSHIFT",2]*x[9]^df["MELAST",2]    #M_DEMA
  
  return(c(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10))
}

for (i in 1:nrow(dftemp2)){
  
  print(i)
  # for (i in 4:4){
  df<-data.frame(matrix(nrow=8, ncol=2))
  df[1,1]<-dftemp2[i,"XSHIFT"]
  df[1,2]<-dftemp2[i,"ROW_XSHIFT"]
  df[2,1]<-dftemp2[i,"XELAST"]
  df[2,2]<-dftemp2[i,"XELAST"]
  #df[3,1]<-dftemp2[i,"TWA"]
  df[3,1]<-dftemp2[i,"Tnew"]
  #df[3,2]<-dftemp2[i,"W_TWA"]
  df[3,2]<-dftemp2[i,"ROW_Tnew"]
  df[4,1]<-dftemp2[i,"Ashift"]
  df[4,2]<-dftemp2[i,"Ashift"]
  df[5,1]<-dftemp2[i,"Mshare"]
  df[5,2]<-dftemp2[i,"ROW_Mshare"]
  
  df[6,1]<-dftemp2[i,"MSUBPAR"]
  df[6,2]<-dftemp2[i,"MSUBPAR"]
  
  df[7,1]<-dftemp2[i,"MSHIFT"]
  df[7,2]<-dftemp2[i,"MSHIFT"]
  
  df[8,1]<-dftemp2[i,"MELAST"]
  df[8,2]<-dftemp2[i,"MELAST"]
  
  rownames(df)<-c("XSHIFT","XELAST","T","ASHIFT","MSHARE","MSUBPAR","MSHIFT","MELAST")
  names(df)<-c("Partner","ROW")
  
  
  
  xstart<-c(1,1,1,1,1,1,1,1,1,1)
  
  # mxit <- 10^8*dftemp2[1,"diviser"]
  z <- nleqslv(xstart, fn3, control=list(maxit=10^8,allowSingular=TRUE)
  )
  
  
  results[i,"HS"]<-dftemp2[i,"HS"]
  results[i,"X"]<-z$x[[1]]
  results[i,"PX"]<-z$x[[2]]
  results[i,"ROW_X"]<-z$x[[3]]
  results[i,"ROW_PX"]<-z$x[[4]]
  
  results[i,"M"]<-z$x[[5]]
  results[i,"ROW_M"]<-z$x[[6]]
  results[i,"PM"]<-z$x[[7]]
  results[i,"ROW_PM"]<-z$x[[8]]
  results[i,"PMA1"]<-z$x[[9]]
  
}


df<-full_join(dftemp2,results,by=c("HS"))


df$change_XS<-(df$XSHIFT/(df$XELAST+1)*(df$PX^(df$XELAST+1)-df$PX0^(df$XELAST+1))) * df$diviser #CH_XGFT(P) CIS
df$ROW_change_XS<-(df$ROW_XSHIFT/(df$XELAST+1)*(df$ROW_PX^(df$XELAST+1)-df$ROW_PX0^(df$XELAST+1))) * df$diviser #CH_XGFT(P) ROW


df$MS_change<- df$diviser * (ifelse((df$MELAST+1) != 0, df$MSHIFT/(df$MELAST+1), 0)*(df$PMA^(df$MELAST+1)-df$PMA1^(df$MELAST+1)))   #CH_MGFT #Welfare effect


tnew <- df1$Tnew[[1]]

df$MELAST <- df$MELAST





#Initial Imports
df$imports <- df$imports *  df$diviser 

df$imp_world <- df$imp_world *  df$diviser  #current import ROW

#New imports
df$M <- df$M * df$diviser #new import 

df$ROW_M <- df$ROW_M * df$diviser #new import ROW

df$change_import <- df$M - df$imports #change in import from country 

df$ROW_change_import <- df$ROW_M - df$imp_world # change in import ROW

#tariff revenue
df$TR_change <- ((df$imports * (df$Tnew/100 - df$TWA/100 )) + 
                   (df$ROW_change_import * df$ROW_TWA/100)) 


df$TR0 <- df$TWA/100 * (df$imports + df$imp_world)

df$TR1 <-  df$TR0 + df$TR_change