Conjugate Transpose - an overview | ScienceDirect Topics

Here is the implementation in python of the algorithm in this article:

#! /usbin/python #---------- # This unusual and intriguing algorithm was originally invented # by Michael V. Klibanov, Professor, Department of Mathematics and Statistics, # University of North Carolina at Charlotte. It is published in the following # paper: # M.V. Klibanov, A.V. Kuzhuget and K.V. Golubnichiy, # "An ill-posed problem for the Black-Scholes equation # for a profitable forecast of prices of stock options on real market data", # Inverse Problems, 32 (2016) 015010. #---------- # Script assumes it's called by crontab, at the opening of the market #----- import numpy as np import pause, datetime from bs4 import BeautifulSoup import requests # Quadratic interpolation of the bid and ask option prices, and linear interpolation in between (https://people.math.sc.edu/kellerlv/Quadratic_Interpolation.pdf) def funcQuadraticInterpolationCoef(values): # There is 'scipy.interpolate.interp1d' too y = np.array(values) A = np.array([[1,0,0],[1,-1,1],[1,-2,4]]) return np.linalg.solve(A,y) # https://en.wikipedia.org/wiki/Polynomial_regression def funcUab(t,coef): return coef[2]*t**2 + coef[1]*t + coef[0] def funcF(s, sa, sb, ua, ub): return (s-sb)*(ua-ub)/(sa-sb) + ub # Initialize the volatility and option lists of 3 values optionBid = [0] # dummy value to pop in the loop optionAsk = [0] # dummy value to pop in the loop volatility = [0] # dummy value to pop in the loop # Initalization for the loop Nt = 4 # even number greater than 2: 4, 6, ... Ns = 2 # even number greater than 0: 2, 4, ... twotau = 2 # not a parameter... alpha = 0.01 # not a parameter... dt = twotau / Nt # time grid step dimA = ( (Nt+1)*(Ns+1), (Nt+1)*(Ns+1) ) # Matrix A dimensions dimb = ( (Nt+1)*(Ns+1), 1 ) # Vector b dimensions A = np.zeros( dimA ) # Matrix A b = np.zeros( dimb ) # Vector b portfolio = 1000000 # Money 'available' securityMargin = 0.00083 # EMPIRICAL: needs to be adjusted when taking into account the transaction fees (should rise, see the article p.8) # Wait 10mn after the opening of the market datet = datetime.datetime.now() datet = datetime.datetime(datet.year, datet.month, datet.day, datet.hour, datet.minute + 10) pause.until(datet) # Record the stock and option values and wait 10mn more def funcRetrieveStockOptionVolatility(): # Stock stock_data_url = "https://finance.yahoo.com/quote/MSFT?p=MSFT" stock_data_html = requests.get(data_url).content stock_content = BeautifulSoup(stock_data_html, "html.parser") stock_bid = content.find("td", {'class': 'Ta(end) Fw(600) Lh(14px)', 'data-test': "BID-value"}) print(stock_bid) stock_ask = content.find("td", {'class': 'Ta(end) Fw(600) Lh(14px)', 'data-test': "ASK-value"}) print(stock_ask) stockOptVol[0] = stock_bid.text.split()[0] stockOptVol[1] = stock_ask.text.split()[0] # Option option_data_url = "https://finance.yahoo.com/quote/MSFT/options?p=MSFT&date=1631836800" option_data_html = requests.get(option_data_url).content option_content = BeautifulSoup(option_data_html, "html.parser") call_option_table = content.find("table", {'class': 'calls W(100%) Pos(r) Bd(0) Pt(0) list-options'}) calls = call_option_table.find_all("tr")[1:] it = 0 for call_option in calls: it+=1 print("it = ", it) if "in-the-money " in str(call_option): itm_calls.append(call_option) print("in the money") itm_put_data = [] for td in BeautifulSoup(str(itm_calls[-1]), "html.parser").find_all("td"): itm_put_data.append(td.text) print(itm_put_data) if itm_put_data[0] == 'MSFT210917C00220000': # One single option stockOptVol[2] = float(itm_put_data[4]) stockOptVol[3] = float(itm_put_data[5]) stockOptVol[4] = float(itm_put_data[-1].strip('%')) else: otm_calls.append(call_option) print("out the money") print("bid = ", option_bid, "\nask = ", option_ask, "\nvol = ",option_vol) return stockOptVol # Record option and volatility stockOptVol = funcRetrieveStockOptionVolatility() optionBid.append(stockOptVol[2]) optionAsk.append(stockOptVol[3]) optionVol.append(stockOptVol[4]) # Wait another 10mn to record a second value for the quadratic interpolation datet = datetime.datetime.now() datet = datetime.datetime(datet.year, datet.month, datet.day, datet.hour, datet.minute + 10) pause.until(datet) stockOptVol = funcRetrieveStockOptionVolatility() optionBid.append(stockOptVol[2]) optionAsk.append(stockOptVol[3]) optionVol.append(stockOptVol[4]) tradeAtTimeTau = False tradeAtTimeTwoTau = False # Run the loop until 30mn before closure datet = datetime.datetime.now() datetend = datetime.datetime(datet.year, datet.month, datet.day, datet.hour + 6, datet.minute + 10) while datet <= datetend: datet = datetime.datetime(datet.year, datet.month, datet.day, datet.hour, datet.minute + 10) optionBid.pop(0) optionAsk.pop(0) optionVol.pop(0) stockOptVol = funcRetrieveStockOptionVolatility() stockBid = stockOptVol[0] stockAsk = stockOptVol[1] optionBid.append(stockOptVol[2]) optionAsk.append(stockOptVol[3]) optionVol.append(stockOptVol[5]) # Trade if required if tradeAtTimeTau == True or tradeAtTimeTwoTau == True: # sell if tradeAtTimeTau == True: portfolio += min(optionAsk[2],sellingPriceAtTimeTau) * 140 # sell 140 options bought 10mn ago tradeAtTimeTau = tradeAtTimeTwoTau sellingPriceAtTimeTau = sellingPriceAtTimeTwoTau sellingPriceAtTimeTwoTau = false else: # forecast the option when no trading # Interpolation coefa = funcQuadraticInterpolationCoef(optionAsk) # quadratic interpolation of the option ask price coefb = funcQuadraticInterpolationCoef(optionBid) # quadratic interpolation of the option bid price coefs = funcQuadraticInterpolationCoef(optionVol) # quadratic interpolation of the volatility sigma sa = stockAsk # stock ask price sb = stockBid # stock bid price ds = (sa - sb) / Ns # stock grid step for k in range (0, Ns+1): # fill the matrix and the vector for j in range (0, Nt+1): Atemp = np.zeros( dimA ) btemp = np.zeros( dimb ) print("k = {k}, j = {j}".format(k=k,j=j)) if k == 0: Atemp[ k*(Nt+1)+j, k*(Nt+1)+j ] = 1 btemp[ k*(Nt+1)+j ] = funcUab(j*dt,coefb) elif k == Ns: Atemp[ k*(Nt+1)+j, k*(Nt+1)+j ] = 1 btemp[ k*(Nt+1)+j ] = funcUab(j*dt,coefa) elif j == 0: Atemp[ k*(Nt+1)+j, k*(Nt+1)+j ] = 1 btemp[ k*(Nt+1)+j ] = funcF( k*ds+sb, sa, sb, funcUab(j*dt,coefa), funcUab(j*dt,coefb) ) elif j == Nt: # do nothing pass else: # main case akj = 0.5*(255*13*3)* funcUab(j*dt, coefs)**2 * (k*ds + sb)**2 dts = (twotau-dt)/Nt * (sa-sb-ds)/Ns #---------- #----- Integral of the generator L #---------- #----- time derivative #---------- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j+1) ] = dts / dt**2 # k,j+1 ~ k,j+1 Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j-1) ] = dts / dt**2 # k,j-1 ~ k,j-1 #----- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j-1) ] = - dts / dt**2 # k,j+1 ~ k,j-1 Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j+1) ] = - dts / dt**2 # k,j-1 ~ k,j+1 #---------- #----- stock derivative #---------- Atemp[ (k+1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] = akj**2 * dts / ds**4 # k+1,j ~ k+1,j Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] = 4 * akj**2 * dts / ds**4 # k,j ~ k,j Atemp[ (k-1)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] = akj**2 * dts / ds**4 # k-1,j ~ k-1,j #----- Atemp[ (k+1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] = -2 * akj**2 * dts / ds**4 # k+1,j ~ k,j Atemp[ (k+0)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] = -2 * akj**2 * dts / ds**4 # k,j ~ k+1,j #----- Atemp[ (k-1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] = -2 * akj**2 * dts / ds**4 # k-1,j ~ k,j Atemp[ (k+0)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] = -2 * akj**2 * dts / ds**4 # k,j ~ k-1,j #----- Atemp[ (k+1)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] = akj**2 * dts / ds**4 # k+1,j ~ k-1,j Atemp[ (k-1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] = akj**2 * dts / ds**4 # k-1,j ~ k+1,j #---------- #----- time and stock derivatives #---------- Atemp[ (k+0)*(Nt+1)+(j+1), (k+1)*(Nt+1)+(j+0) ] = akj * dts / (dt*ds**2) # k,j+1 ~ k+1,j Atemp[ (k+1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+1) ] = akj * dts / (dt*ds**2) # k+1,j ~ k,j+1 #----- Atemp[ (k+0)*(Nt+1)+(j-1), (k+1)*(Nt+1)+(j+0) ] = - akj * dts / (dt*ds**2) # k,j-1 ~ k+1,j Atemp[ (k+1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j-1) ] = - akj * dts / (dt*ds**2) # k+1,j ~ k,j-1 #---------- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j+0) ] = -2 * akj * dts / (dt*ds**2) # k,j+1 ~ k,j Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+1) ] = -2 * akj * dts / (dt*ds**2) # k,j ~ k,j+1 #----- Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j+0) ] = 2 * akj * dts / (dt*ds**2) # k,j-1 ~ k,j Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j-1) ] = 2 * akj * dts / (dt*ds**2) # k,j ~ k,j-1 #---------- Atemp[ (k+0)*(Nt+1)+(j+1), (k-1)*(Nt+1)+(j+0) ] = akj * dts / (dt*ds**2) # k,j+1 ~ k-1,j Atemp[ (k-1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+1) ] = akj * dts / (dt*ds**2) # k-1,j ~ k,j+1 #----- Atemp[ (k+0)*(Nt+1)+(j-1), (k-1)*(Nt+1)+(j+0) ] = - akj * dts / (dt*ds**2) # k,j-1 ~ k-1,j Atemp[ (k-1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j-1) ] = - akj * dts / (dt*ds**2) # k-1,j ~ k,j-1 #---------- #---------- #----- Regularisation term - using alpha = 0.01 #---------- #---------- #----- H2 norm: 0 derivative #---------- Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] += alpha # k,j ~ k,j #----- coef = funcF( k*ds+sb, sa, sb, funcUab(j*dt,coefa), funcUab(j*dt,coefb) ) btemp[ (k+0)*(Nt+1)+(j+0) ] += alpha * 2 * coef #---------- #----- H2 norm: time derivative #---------- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j+1) ] += alpha / dt**2 # k,j+1 ~ k,j+1 Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j-1) ] += alpha / dt**2 # k,j-1 ~ k,j-1 #----- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j-1) ] += -alpha / dt**2 # k,j+1 ~ k,j-1 Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j+1) ] += -alpha / dt**2 # k,j-1 ~ k,j+1 #----- coef = ( funcF( k*ds+sb, sa, sb, funcUab((j+1)*dt,coefa), funcUab((j+1)*dt,coefb) ) \ - funcF( k*ds+sb, sa, sb, funcUab((j-1)*dt,coefa), funcUab((j-1)*dt,coefb) ) ) / dt btemp[ (k+0)*(Nt+1)+(j+1) ] += alpha * 2 * coef btemp[ (k+0)*(Nt+1)+(j-1) ] += - alpha * 2 * coef #---------- #----- H2 norm: stock derivative #---------- Atemp[ (k+1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] += alpha / ds**2 # k+1,j ~ k+1,j Atemp[ (k-1)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] += alpha / ds**2 # k-1,j ~ k-1,j #----- Atemp[ (k+1)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] += -alpha / ds**2 # k+1,j ~ k-1,j Atemp[ (k-1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] += -alpha / ds**2 # k-1,j ~ k+1,j #----- coef = ( funcUab(j*dt,coefa) - funcUab(j*dt,coefb) ) / (sa - sb) btemp[ (k+1)*(Nt+1)+(j+0) ] += alpha * 2 * coef btemp[ (k-1)*(Nt+1)+(j+0) ] += - alpha * 2 * coef #---------- #----- H2 norm: stock and time derivative #---------- Atemp[ (k+1)*(Nt+1)+(j+1), (k+1)*(Nt+1)+(j+1) ] += alpha / (ds*dt) # k+1,j+1 ~ k+1,j+1 Atemp[ (k-1)*(Nt+1)+(j+1), (k-1)*(Nt+1)+(j+1) ] += alpha / (ds*dt) # k-1,j+1 ~ k-1,j+1 Atemp[ (k-1)*(Nt+1)+(j-1), (k-1)*(Nt+1)+(j-1) ] += alpha / (ds*dt) # k-1,j-1 ~ k-1,j-1 Atemp[ (k+1)*(Nt+1)+(j-1), (k+1)*(Nt+1)+(j-1) ] += alpha / (ds*dt) # k+1,j-1 ~ k+1,j-1 #---------- Atemp[ (k+1)*(Nt+1)+(j+1), (k-1)*(Nt+1)+(j+1) ] += -alpha / (ds*dt) # k+1,j+1 ~ k-1,j+1 Atemp[ (k+1)*(Nt+1)+(j+1), (k+1)*(Nt+1)+(j-1) ] += -alpha / (ds*dt) # k+1,j+1 ~ k+1,j-1 Atemp[ (k+1)*(Nt+1)+(j+1), (k-1)*(Nt+1)+(j-1) ] += alpha / (ds*dt) # k+1,j+1 ~ k-1,j-1 #----- Atemp[ (k-1)*(Nt+1)+(j+1), (k+1)*(Nt+1)+(j+1) ] += -alpha / (ds*dt) # k-1,j+1 ~ k+1,j+1 Atemp[ (k+1)*(Nt+1)+(j-1), (k+1)*(Nt+1)+(j+1) ] += -alpha / (ds*dt) # k+1,j-1 ~ k+1,j+1 Atemp[ (k-1)*(Nt+1)+(j-1), (k+1)*(Nt+1)+(j+1) ] += alpha / (ds*dt) # k-1,j-1 ~ k+1,j+1 #---------- Atemp[ (k-1)*(Nt+1)+(j+1), (k+1)*(Nt+1)+(j-1) ] += alpha / (ds*dt) # k-1,j+1 ~ k+1,j-1 Atemp[ (k-1)*(Nt+1)+(j+1), (k-1)*(Nt+1)+(j-1) ] += -alpha / (ds*dt) # k-1,j+1 ~ k-1,j-1 #----- Atemp[ (k+1)*(Nt+1)+(j-1), (k-1)*(Nt+1)+(j+1) ] += alpha / (ds*dt) # k+1,j-1 ~ k-1,j+1 Atemp[ (k-1)*(Nt+1)+(j-1), (k-1)*(Nt+1)+(j+1) ] += -alpha / (ds*dt) # k-1,j-1 ~ k-1,j+1 #---------- Atemp[ (k+1)*(Nt+1)+(j-1), (k-1)*(Nt+1)+(j-1) ] += -alpha / (ds*dt) # k+1,j-1 ~ k-1,j-1 #----- Atemp[ (k-1)*(Nt+1)+(j-1), (k+1)*(Nt+1)+(j-1) ] += -alpha / (ds*dt) # k-1,j-1 ~ k+1,j-1 #---------- coef = ( funcUab((j+1)*dt,coefa) - funcUab((j+1)*dt,coefb) \ - funcUab((j-1)*dt,coefa) + funcUab((j-1)*dt,coefb) ) / (dt * (sa - sb)) btemp[ (k+1)*(Nt+1)+(j+1) ] += alpha * 2 * coef / (ds*dt) btemp[ (k-1)*(Nt+1)+(j+1) ] += - alpha * 2 * coef / (ds*dt) btemp[ (k-1)*(Nt+1)+(j-1) ] += - alpha * 2 * coef / (ds*dt) btemp[ (k+1)*(Nt+1)+(j-1) ] += alpha * 2 * coef / (ds*dt) #---------- #----- H2 norm: stock second derivative #---------- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j+1) ] += alpha / dt**4 # k,j+1 ~ k,j+1 Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] += 4 * alpha / dt**4 # k,j ~ k,j Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j-1) ] += alpha / dt**4 # k,j-1 ~ k,j-1 #----- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j+0) ] += -2 * alpha / dt**4 # k,j+1 ~ k,j Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+1) ] += -2 * alpha / dt**4 # k,j ~ k,j+1 #----- Atemp[ (k+0)*(Nt+1)+(j+1), (k+0)*(Nt+1)+(j-1) ] += alpha / dt**4 # k,j+1 ~ k,j-1 Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j+1) ] += alpha / dt**4 # k,j-1 ~ k,j+1 #----- Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j-1) ] += -2 * alpha / dt**4 # k,j ~ k,j-1 Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j+0) ] += -2 * alpha / dt**4 # k,j-1 ~ k,j #---------- #----- H2 norm: time second derivative #---------- Atemp[ (k+1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] += alpha / ds**4 # k+1,j ~ k+1,j Atemp[ (k+0)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] += 4 * alpha / ds**4 # k,j ~ k,j Atemp[ (k+1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] += alpha / ds**4 # k-1,j ~ k-1,j #----- Atemp[ (k+1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] += -2 * alpha / ds**4 # k+1,j ~ k,j Atemp[ (k+0)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] += -2 * alpha / ds**4 # k,j ~ k+1,j #----- Atemp[ (k+1)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] += alpha / ds**4 # k,j ~ k,j Atemp[ (k-1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] += alpha / ds**4 # k,j ~ k,j #----- Atemp[ (k+0)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] += -2 * alpha / ds**4 # k,j ~ k-1,j Atemp[ (k-1)*(Nt+1)+(j+0), (k+0)*(Nt+1)+(j+0) ] += -2 * alpha / ds**4 # k-1,j ~ k,j #---------- coef = ( funcF( k*ds+sb, sa, sb, funcUab((j+1)*dt,coefa), funcUab((j+1)*dt,coefb) ) \ - 2 * funcF( k*ds+sb, sa, sb, funcUab((j+0)*dt,coefa), funcUab((j+0)*dt,coefb) ) \ + funcF( k*ds+sb, sa, sb, funcUab((j-1)*dt,coefa), funcUab((j-1)*dt,coefb) ) ) / dt**2 btemp[ (k+0)*(Nt+1)+(j+1) ] += alpha * 2 * coef / dt**2 btemp[ (k+0)*(Nt+1)+(j+0) ] += - alpha * 4 * coef / dt**2 btemp[ (k+0)*(Nt+1)+(j-1) ] += alpha * 2 * coef / dt**2 #---------- #---------- #----- Boundary de-computation #---------- if k+1 == Ns: Atemp[ (k+1)*(Nt+1)+(j+0), (k+1)*(Nt+1)+(j+0) ] = 0 # k+1,j ~ k+1,j Atemp[ (k+1)*(Nt+1)+(j+1), (k+1)*(Nt+1)+(j+1) ] = 0 # k+1,j+1 ~ k+1,j+1 Atemp[ (k+1)*(Nt+1)+(j-1), (k+1)*(Nt+1)+(j-1) ] = 0 # k+1,j-1 ~ k+1,j-1 btemp[ (k+1)*(Nt+1)+(j+0) ] = 0 # k+1,j btemp[ (k+1)*(Nt+1)+(j+1) ] = 0 # k+1,j+1 btemp[ (k+1)*(Nt+1)+(j-1) ] = 0 # k+1,j-1 if k-1 == 0: Atemp[ (k-1)*(Nt+1)+(j+0), (k-1)*(Nt+1)+(j+0) ] = 0 # k-1,j ~ k-1,j Atemp[ (k-1)*(Nt+1)+(j+1), (k-1)*(Nt+1)+(j+1) ] = 0 # k-1,j+1 ~ k-1,j+1 Atemp[ (k-1)*(Nt+1)+(j-1), (k-1)*(Nt+1)+(j-1) ] = 0 # k-1,j-1 ~ k-1,j-1 btemp[ (k-1)*(Nt+1)+(j+0) ] = 0 # k-1,j btemp[ (k-1)*(Nt+1)+(j+1) ] = 0 # k-1,j+1 btemp[ (k-1)*(Nt+1)+(j-1) ] = 0 # k-1,j-1 if j-1 == 0: Atemp[ (k+0)*(Nt+1)+(j-1), (k+0)*(Nt+1)+(j-1) ] = 0 # k,j-1 ~ k,j-1 Atemp[ (k+1)*(Nt+1)+(j-1), (k+1)*(Nt+1)+(j-1) ] = 0 # k+1,j-1 ~ k+1,j-1 Atemp[ (k-1)*(Nt+1)+(j-1), (k-1)*(Nt+1)+(j-1) ] = 0 # k-1,j-1 ~ k-1,j-1 btemp[ (k+0)*(Nt+1)+(j-1) ] = 0 # k,j-1 btemp[ (k+1)*(Nt+1)+(j-1) ] = 0 # k+1,j-1 btemp[ (k-1)*(Nt+1)+(j-1) ] = 0 # k-1,j-1 #---------- pass print("-----") print("Atemp = ") print(Atemp) print("-----") print("btemp = ") print(btemp) print("-----") print("-----") A = A + Atemp b = b + btemp print("-----") print("A = ") print(A) print("-----") print("b = ") print(b) print("-----") print("-----") input("Press Enter to continue...") # Conjugate gradient algorithm: https://en.wikipedia.org/wiki/Conjugate_gradient_method x = np.zeros(N).reshape(N,1) r = b - np.matmul(A,x) p = r rsold = np.dot(r.transpose(),r) for i in range(len(b)): Ap = np.matmul(A,p) alpha = rsold / np.matmul(p.transpose(),Ap) x = x + alpha * p r = r - alpha * Ap rsnew = np.dot(r.transpose(),r) if np.sqrt(rsnew) < 1e-16: break p = r + (rsnew / rsold) * p rsold = rsnew print("it = ", i) print("rsold = ", rsold) # Trading strategy sm = (sa + sb)/2 if x[Ns/2*(Nt+1)+Nt/2] >= optionAsk[0] + securityMargin: tradeAtTimeTau = True sellingPriceAtTimeTau = x[Ns/2*(Nt+1)+Nt/2] portfolio -= 140 * optionAsk # buy 140 options if x[Ns/2*(Nt+1)+Nt] >= optionAsk[0] + securityMargin: tradeAtTimeTwoTau = True sellingPriceAtTimeTwoTau = x[Ns/2*(Nt+1)+Nt] portfolio -= 140 * optionAsk # buy 140 options pause.until(datet) # Wait 10mn before the next loop pause.until(datet) datet = datetime.datetime.now() # Time should be around 20mn before closure datet = datetime.datetime(datet.year, datet.month, datet.day, datet.hour, datet.minute + 10) if tradeAtTimeTau == True: # sell stockOptVol = funcRetrieveStockOptionVolatility() optionAsk.pop(0) optionAsk.append(stockOptVol[3]) portfolio += min(optionAsk[2],sellingPriceAtTimeTau) * 140 # Wait 10mn more to sell the last options pause.until(datet) # it should be around 10mn before closure if tradeAtTimeTwoTau == True: # sell stockOptVol = funcRetrieveStockOptionVolatility() optionAsk.pop(0) optionAsk.append(stockOptVol[3]) portfolio += min(optionAsk[2],sellingPriceAtTimeTwoTau) * 140 # Market closure 

Don't put money on this as I'm still debugging (I bet you half a bitcoin I have mistaken a few indices in the H_2 norm)... Here is the discretisation formula I used, to copy-paste on latexbase:

\documentclass[12pt]{article} \usepackage{amsmath} \usepackage[latin1]{inputenc} \title{Klibanov algorithm} \author{Discretisation formula} \date{\today} \begin{document} \maketitle Let $$ a_{k,j} = \frac12\sigma(j\delta_\tau)^2\times(255\times13\times3)\times(k\delta_s+s_a)^2, $$ then \begin{alignat*}{3} J_\alpha(u) = & \sum_{k=1}^{N_s} \sum_{j=1}^{N_t} \left| \frac{u_{k,j+1} - u_{k,j-1}}{\delta_\tau} + a_{k,j} \frac{u_{k+1,j} - 2u_{k,j} + u_{k-1,j}}{\delta_s^2}\right|^2\frac{2\tau - \delta_\tau}{N_t}\frac{s_a - s_b - \delta_s}{N_s}\\ & + \alpha \sum_{k=1}^{N_s} \sum_{j=1}^{N_t} \left| u_{k,j} - F_{k,j}\right|^2 \\ & \qquad + \left| \frac{u_{k,j+1} - u_{k,j-1}}{\delta_t} - \frac{F_{k,j+1} - F_{k,j-1}}{\delta_t}\right|^2 \\ & \qquad + \left| \frac{u_{k+1,j} - u_{k-1,j}}{\delta_s} - \frac{u_{a,j} - u_{b,j}}{s_a - s_b}\right|^2 \\ & \qquad + \left| \frac{(u_{k+1,j+1} - u_{k-1,j+1}) - (u_{k+1,j-1} - u_{k-1,j-1})}{\delta_s\delta_t} \right. \\ & \qquad \qquad \left. - \frac{(u_{a,j+1} - u_{b,j+1}) - (u_{a,j-1} - u_{b,j-1})}{(s_a-s_b)\delta_t}\right|^2 \\ & \qquad + \left| \frac{u_{k,j+1} - 2u_{k,j} + u_{k,j-1}}{\delta_\tau^2} - \frac{F_{k,j+1} - 2F_{k,j} + F_{k,j-1}}{\delta_\tau^2} \right|^2 \\ & \qquad + \left| \frac{u_{k+1,j} - 2u_{k,j} + u_{k-1,j}}{\delta_s^2}\right|^2 \end{alignat*} %% \left| \right|^2 with $\tau = 1$ unit of time (for example 10mn). \end{document} 

Let me know if you see something wrong... And if you want to contribute, feel free

Hey y'all. I'm reading Elementary Matrix Algebra by Franz Hohn (from 1962). He insists that two complex vectors X and Y are orthogonal if and only if

X*Y = 0

Where X* is the transposed conjugate of X.
My question I suppose has little to do with this text, but it did lead to the question nonetheless: what is orthogonality in Cⁿ, the n-dimensional complex space (is that even correct notation for it)?
My initial instinct is to use a simple example (using semicolons to end rows, commas to end columns): X = [i; 0; 0]; Y = [1; 0; 0]. Obviously (right?) 1 and i are orthogonal in the complex plane, but X*Y = [-i, 0, 0]*[1;0;0] = (-i + 0 + 0) = -i ≠ 0.
Can someone please explain orthogonality with complex vectors?

Good day,
I am working with large matrices, lets call one A, (generally 500x500, but I dont want to assume the matrix is square) where the elements in the matrix A range from (inslusive) 0 to 255 (they are pixel brightness components, so are also integers). I am applying the SVD to A. I would like to be able to say (and prove, or reference) what number the resulting singular values (which are the square roots of the eigenvalues of AA* (where A* is the conjugate transpose of A)) are definitely less than.
i.e. I want to find b, where;

max{ singular values of A } < b 

I've tried to understand the min-max theorem but have not got far - if this theorem is relevant to me then please say and if possible give an example.
I 'guessed' an assumpiton, that the greatest singular value would come from a matrix containing 255 for every element, and so far I've been unable to prove this assumption wrong - but I've not proven it. I've also shown (again, not proven) that when a matrix (any dimension) contains the same value for every element that the max singular value increases linearly as the matrix dimension increases. This means that the inequality I'm looking for will need to take into account matrix height and width (or perhaps just the matrix rank). I can assume that the max value of the height and width are,

max{width} = max{height} = (2^16)-1 = 65535 

I'm not looking just for the best upper bound (i.e. the lowest upper bound). Any upper bound might be enough.
Any help is obviously appreciated.

Good day,
I am working with large matrices, lets call one A, (generally 500x500, but I dont want to assume the matrix is square) where the elements in the matrix A range from (inslusive) 0 to 255 (they are pixel brightness components, so are also integers). I am applying the SVD to A. I would like to be able to say (and prove, or reference) what number the resulting singular values (which are the square roots of the eigenvalues of AA* (where A* is the conjugate transpose of A)) are definitely less than.
i.e. I want to find b, where;

max{ singular values of A } < b 

I've tried to understand the min-max theorem but have not got far - if this theorem is relevant to me then please say and if possible give an example.
I 'guessed' an assumpiton, that the greatest singular value would come from a matrix containing 255 for every element, and so far I've been unable to prove this assumption wrong - but I've not proven it. I've also shown (again, not proven) that when a matrix (any dimension) contains the same value for every element that the max singular value increases linearly as the matrix dimension increases. This means that the inequality I'm looking for will need to take into account matrix height and width (or perhaps just the matrix rank). I can assume that the max value of the height and width are,

max{width} = max{height} = (2^16)-1 = 65535 

I'm not looking just for the best upper bound (i.e. the lowest upper bound). Any upper bound might be enough.
Any help is obviously appreciated.

I heard that linear algebra especially matrix algebra including singular value decomposition, symmetric, hermitian, conjugate transpose, unitary geometry, transposes, and spectral theory show up in multivariate statistics and stochastic processes. Is multivariate statistics really dependent on understanding advanced linear algebra? Should I avoid multivariate statistics and stochastic processes if the advanced linear algebra and diagonalization methods of matricies don't make sense to me?
Whenever I read a stochastic processes book or book on multivariate statistics, the matrix algebra used when doing change of basis scares me and confuses me.
I've also look at a book on probability theory which is a prerequisite for stochastic processes which is applied calculus and all the computations or proofs on for example leverage and cook's distance seem messy and not elegant. I don't think I'd like mathematical stats.

I know they generalize the conjugate transpose of a matrix but I am not sure how. I also know they involve a restriction mapping. I have seen some examples of them but I don't understand why the examples satisfy the definition.

FYI: VS2013, C++, using glm's gtx fdualquat but I defined my own mat4 cast because I couldn't find a glm::mat4_cast(...) that was overloaded for a glm::fdualquat. I am changing all the signs on my controls for the moment so that I can keep working on it, but I'd appreciate any suggestions on why my axes all flipped around (note: this also had the effect of reversing all my rotations).
Attempts to fix the problem:

• I tried calculating the transpose of the matrix in the event that I mixed up row-first and column-first convention, but that that destroyed my model_to_world matrix :), so that's not the problem.
  
• I also tried changing from left-multiplying the fdualquat's to right-multiplying them, but that didn't change the reversed axis problem.
  

I got the code conversion from this paper: "A Beginners Guide to Dual-Quaternions" (http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.407.9047, but the download link at this site seems down at the moment). The example code appeared to be in C#, but I didn't see a mention of the language, so I don't know for sure. I have double and triple checked my code with the reference and I'm sure it's correct.
Here's the code:

glm::mat4 dual_quat_to_mat4(const glm::fdualquat &dq) { glm::fdualquat norm_dq = glm::normalize(dq); // get an identity matrix glm::mat4 mat; // extract rotational information float scalar = norm_dq.real.w; float x = norm_dq.real.x; float y = norm_dq.real.y; float z = norm_dq.real.z; mat[0][0] = (scalar * scalar) + (x * x) - (y * y) - (z * z); mat[0][1] = (2 * x * y) + (2 * scalar * z); mat[0][2] = (2 * x * z) - (2 * scalar * y); mat[1][0] = (2 * x * y) - (2 * scalar * z); mat[1][1] = (scalar * scalar) + (y * y) - (x * x) - (z * z); mat[1][2] = (2 * y * z) + (2 * scalar * x); mat[2][0] = (2 * x * z) + (2 * scalar * y); mat[2][1] = (2 * y * z) - (2 * scalar * x); mat[2][2] = (scalar * scalar) + (z * z) - (x * x) - (y * y); // extract translational information glm::fquat trans = (norm_dq.dual * 2.0f) * glm::conjugate(norm_dq.real); mat[3][0] = trans.x; mat[3][1] = trans.y; mat[3][2] = trans.z; // the last row remains untouched; that row is the realm of clip space and perspective division return mat; } 

Again, my problem is that all my self-defined standard vectors (defined to establish a consistent set of unit vectors) flipped around. My world up vector of <0.0f, +1.0f, 0.0f> became <0.0f, -1.0f, 0.0f>, my world left vector of <-1.0f, 0.0f, 0.0f> became <+1.0f, 0.0f, 0.0f>, and my world forward vector of <0.0f, 0.0f, -1.0f> became <0.0f, 0.0f, +1.0f>.
Ex: I try setting an object to be at the origin, but +8.0f up in the air. That translation used to move it up above my ground plane, but after I switched to dual quats +8.0f moves it down 8.0f below the plane.